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ABSTRACT 

We simulate the assembly of a massive rich cluster and the formation of its constituent 
galaxies in a flat, low-density universe. Our most accurate model follows the collapse, 
the star-formation history and the orbital motion of all galaxies more luminous than 
the Fornax dwarf spheroidal, while dark halo structure is tracked consistently through- 
out the cluster for all galaxies more luminous than the SMC. Within its virial radius 
this model contains about 2 x 10 7 dark matter particles and almost 5000 distinct 
dynamically resolved galaxies. Simulations of this same cluster at a variety of resolu- 
tions allow us to check explicitly for numerical convergence both of the dark matter 
structures produced by our new parallel N-body and substructure identification codes, 
and of the galaxy populations produced by the phenomenological models we use to 
follow cooling, star formation, feedback and stellar aging. This baryonic modelling is 
tuned so that our simulations reproduce the observed properties of isolated spirals 
outside clusters. Without further parameter adjustment our simulations then produce 
a luminosity function, a mass-to-light ratio, luminosity, number and velocity disper- 
sion profiles, and a morphology-radius relation which are similar to those observed in 
real clusters. In particular, since our simulations follow galaxy merging explicitly, we 
can demonstrate that it accounts quantitatively for the observed cluster population of 
bulges and elliptical galaxies. 

Key words: galaxies: formation - galaxies: clusters: general - dark matter. 
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1 INTRODUCTION 

The last two decades have witnessed substantial progress 
towards an understanding of hierarchical galaxy formation 
within the framework of a universe dominated by cold dark 
matter (CDM). For an appropriate choice of the cosmologi- 
cal parameters, the CDM theory provides a remarkably suc- 
cessful description of large-scale structure formation, and it 
is in good agreement with a large variety of observational 
data. Much of this progress has been achieved by detailed 
analytical and numerical studies of the collisionless dynam- 
ics of the dark matter. As a result, this part of cosmic evo- 
lution is now quite well understood. However, the actual 
formation of the luminous parts of galaxies within CDM uni- 
verses involves many complex physical processes in addition 
to gravity, for example, shocking and cooling of gas, and 
star formation with its attendant regulation and feedback 
mechanisms. The theoretical modeling of important aspects 
of these processes is still highly uncertain. 
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Not surprisingly, the lack of precise specifications for 
treating the relevant physics has also hampered direct nu- 
merical studies of galaxy formation. In addition, such studies 
are confronted with the huge range of scales over which these 
physical processes interact. Difficulties in finding robust and 
appropriate algorithms for handling "subgrid" physics have 
so far prevented hydrodynamical simulations from reproduc- 
ing many basic properties of galaxies, although more recent 
work is beginning to achieve some notable successes both 
for individual galaxies (e.g. Katz & Gunn 1991; Navarro & 
White 1994; Steinmetz & Miiller 1995; Mihos & Hernquist 
1996; Walker et al. 1996; Navarro & Steinmetz 1997; Stein- 
metz & Navarro 1999) and for the distribution of galaxy 
populations (e.g. Cen & Ostriker 1993, 2000; Katz et al. 
1996, 1999; Weinberg et al. 1997, 2000; Steinmetz & Miiller 
1995; Blanton et al. 1999; Pearce et al. 1999, 2000). 

Much of our current understanding of galaxy formation 
has been learned through 'semi-analytic' techniques, as laid 
out originally by White & Frenk (1991), Cole (1991) and 
Lacey & Silk (1991), building on the scenario first sketched 
by White & Rees (1978) . In these models, each of the compli- 
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cated and interacting physical processes involved in galaxy 
formation is approximated using a simplified, physically 
based model. These processes include the growth through 
accretion and merging of dark matter haloes, the shock heat- 
ing and virialization of gas within these haloes, the radiative 
cooling of gas and its settling to a rotationally supported 
disk, star formation, and the resulting feedback from super- 
novae and stellar winds, the evolution of stellar populations, 
absorption and reradiation by dust, and galaxy merging with 
its accompanying starbursts and morphological transforma- 
tions. At the expense of uncertainties introduced by the sim- 
plifying assumptions, semi-analytic techniques can access a 
much larger dynamic range than numerical simulations, they 
allow a fast exploration of parameter space and of the in- 
fluence of the specific simplifying assumptions chosen, and 
they facilitate direct comparisons with a wide range of ob- 
servational data. 

Over the last few years, a number of groups have used 
semi-analytical models to study galaxy formation, and to in- 
terpret observations of galaxy populations at low and high 
redshift (Lacey et al. 1993; Kauffmann et al. 1993, 1994; Cole 
et al. 1994; Heyl et al. 1995; Baugh et al. 1996a,b; Kauff- 
mann 1995a,b, 1996a,b; Guiderdoni et al. 1998; Kauffmann 
& Chariot 1998; Baugh et al. 1998, 1999; Mo et al. 1998, 
1999; Devriendt et al. 1998; Somerville & Primack 1999; 
Kauffmann & Haehnelt 2000; Cole et al. 2000; van den Bosch 
2000; Boisser & Prantzos 2000; Haehnelt & Kauffmann 2000; 
Somerville et al. 2000). Stellar population synthesis mod- 
elling allows detailed photometric comparison with observa- 
tion, including studies of the strong evolution apparent in 
observed high redshift galaxy samples. With a small num- 
ber of free parameters, semi-analytic models have been quite 
successful in allowing a unified and coherent interpretation 
of a broad range of galaxy properties, for example, lumi- 
nosity functions, Tully-Fisher and Faber- Jackson relations, 
number counts, the distributions of morphology, color and 
size, global and individual star formation histories, back- 
ground radiation contributions from the UV to the far IR, 
clustering strengths, and the observed relations between 
AGN's and their host galaxies. Most studies concentrate on a 
few specific issues. Oversimplifications in the adopted physi- 
cal models often show up as inconsistencies with observation 
in other areas. For example, until recently most studies had 
difficulty simultaneously to fit the zero-point of the Tully- 
Fisher relation and the luminosity function of field galaxies. 
Recent work has removed some of the most serious over- 
simplifications in this area and has reduced the discrepancy 
significantly (Somerville & Primack 1999; van den Bosch 
2000). 

The construction of dark matter merging history trees 
(Kauffmann & White 1993; Somerville & Kolatt 1999; Sheth 
& Lemson 1999) is an important ingredient in semi-analytic 
models. In most studies, Monte-Carlo realizations of merg- 
ing histories for individual objects are generated using 
the extended Press-Schechter formalism (Press & Schechter 
1974; Bond et al. 1991; Lacey & Cole 1993). A disadvantage 
of this approach is that there is little information about the 
spatial distribution of galaxies, although two point corre- 
lations can be estimated using the methods introduced by 
Mo & White (1996) (see also Baugh et al. 1999). In or- 
der to study clustering in more detail semi-analytic models 
have been combined with cosmological N-body simulations; 



the galaxy population within each of the virialised haloes 
present in a given output of the simulation is created by a 
Monte Carlo "semi-analytic" realisation of its prior history 
and these galaxies are then attached to the halo centre and 
to random "satellite" particles within the halo (Kauffmann 
et al. 1997; Governato et al. 1998; Benson et al. 2000a,b; 
Wechsler et al. 2000) . This allows mock catalogues of galax- 
ies to be constructed which contain all the spatial and kine- 
matic information of real redshift surveys. 

In a natural extension of this approach, one can use 
N-body simulations not only to provide the mass distribu- 
tion at a given time, but also to reconstruct individual halo 
merging histories, a scheme first tried in a crude form by 
White et al. (1987). This allows one to avoid the uncer- 
tainties inherent in the Press-Schechter formalism, and it 
provides spatial and kinematic distributions not just for the 
final galaxies but for their progenitors at all earlier times 
as well. Thus it is effectively equivalent to a full dynami- 
cal simulation of galaxy formation and clustering, but with 
the advantage that the computationally intensive part of the 
procedure, the original N-body simulation, does not need to 
be repeated every time the assumptions about baryonic pro- 
cesses are changed. The disadvantages in comparison with 
the simpler scheme just discussed are that the resolution of 
the merging history trees is limited by that of the N-body 
simulation and that particle data must be stored sufficiently 
often for the trees to follow adequately the growth of struc- 
ture. The method thus requires frequent data dumps from 
high resolution simulations of cosmologically representative 
volumes, leading to a substantial raw data volume. 

Recently, Roukemaet al. (1997) studied merging history 
trees directly from N-body simulations using scale-free sim- 
ulations and a rather limited number of simulation outputs. 
A much more extensive study has been published by Kauff- 
mann et al. (1999a, hereafter KCDW), who constructed 
merging history trees from two high-resolution N-body sim- 
ulations using a total of 51 output times between redshift 
z = 20 and z = 0. KCDW grafted semi-analytic models of 
galaxy formation onto the simulated merger trees to study 
how the clustering of galaxies is related to intrinsic proper- 
ties like luminosity, colour and morphology. In subsequent 
papers, they used this methodology to predict the evolution 
of clustering to high redshift (Kauffmann et al. 1999b), to 
construct realistically selected mock redshift surveys (Diafe- 
rio et al. 1999), and to study the spatial and kinematic dis- 
tributions of galaxies within clusters (Diaferio et al. 2000). 

Using a somewhat different technique, van Kampen 
et al. (1999) also made use of the full merging history of N- 
body simulations. They modified an existing N-body code 
so that heavier 'tracer' particles were introduced during the 
execution of the simulation. These tracer particles, identi- 
fied with galaxies, replaced locally overdense groups of the 
orginal particles. This approach neglects the internal struc- 
ture of galaxy/halo systems, for example the fact that such 
systems can be stripped of much of their mass as they orbit 
within a cluster. In addition it is somewhat inflexible since 
any change in the assumptions about how baryonic part of 
galaxies forms and evolves requires the simulation to be re- 
peated. 

Our approach follows the methodology of KCDW, ex- 
tending it to deal with higher resolution simulations. In such 
simulations, substantial substructure, partially stripped 
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haloes of cluster galaxies, can be identified within dark mat- 
ter clumps corresponding to galaxy groups and clusters. To 
this end, we study four simulations of the formation of a 
rich cluster of galaxies. These simulations follow the same 
object, a Coma-like cluster of mass 8 x 10 14 /i _1 M0 in a flat, 
low-density universe, but use different mass resolutions. We 
resolve the virial region of the final cluster with 0.12, 0.61, 
3.5 and 20 million particles, respectively, and we sample the 
field in the region immediately around the cluster at the 
same resolution but with about twice as many particles in 
each case. The cosmological tidal field is represented by an 
additional boundary region with ~ 3.1 million particles ex- 
tending to a distance of 70/i _1 Mpc from the cluster. This 
sequence allows us to test explicitly how our results depend 
on numerical resolution. 

We develop a new algorithm, SUBFIND, for identifying 
substructure within the clumps formed in these simulations. 
This algorithm defines 'subhaloes' as locally overdense, self- 
bound particle groups, and is able to detect hierarchies of 
substructure using just a single simulation output. As in 
KCDW, we have stored 51 simulation snapshots from z = 20 
to the current epoch, and we trace the merging history of 
groups and their subhaloes from output to output. We mod- 
ify the semi-analytic recipes employed by KCDW to allow 
the inclusion of subhaloes, and we analyse the changes re- 
sulting from this increase in fidelity to the physical system 
modelled. 

In particular, we study the luminosity function of the 
cluster, its mass-to-light ratio, the Tully-Fisher relation of 
spirals in the surrounding field, the Faber- Jackson relation 
of field and cluster ellipticals, and the B — V colors of our 
model galaxies. We show that the new subhalo-scheme gives 
rise to a pronounced radial segregation by morphology. We 
investigate luminosity segregation by comparing the radial 
profiles of galaxies by number, morphology and luminosity 
with that of the dark matter. Finally, we study how the 
velocity dispersion profile of the galaxies depends on their 
luminosity and colour and compare it with the correspond- 
ing quantity for the dark matter. For all these quantities we 
are able to demonstrate good numerical convergence for all 
galaxies brighter than the SMC. 

Our approach also allows a detailed study of the forma- 
tion history of the cluster and its galaxies. In particular, we 
have analysed the star formation history of field and cluster 
galaxies, the spatial and temporal origin of the 'first' stars 
that end up in the cluster (White & Springel 1999), and the 
evolution of the merger rate of galaxies. These results are 
discussed in a companion paper. 

Interestingly, the new subhalo analysis improves the 
agreement with observational data, especially with respect 
to the cluster luminosity function. Most of the bright galax- 
ies in the final, highly- resolved cluster are still connected 
to well localized subhaloes within the smooth dark matter 
background of the cluster. There is hence no need to esti- 
mate merging timescales within the cluster using dynamical 
friction or approximate cross-section arguments. Mergers are 
treated automatically in a fully dynamically consistent way. 
As we will discuss, inaccuracies in estimated merging times 
can lead to a problem of excessively bright first ranked clus- 
ter galaxies in the simpler methodology which does not track 
subhaloes. This problem goes away in our refined approach 
which also demonstrates explicitly that merging can account 



for the observed fractions of elliptical and bulge-dominated 
galaxies and for their distribution within clusters. 

This paper is organized as follows. In Section 2 we de- 
scribe the N-body simulations, and in Section 3 we review 
the techniques of KCDW, and our specific implementation 
of them. In Section 4 we discuss our techniques for identify- 
ing dark matter substructure within larger haloes, and our 
methods of tracing it from output to ouput in the simula- 
tions. We then describe the implementation of semi-analytic 
models including this subhalo information, and we present 
results obtained with these prescriptions in Section 5. Fi- 
nally, we discuss some aspects of our findings in Section 6. 



2 CLUSTER SIMULATIONS 

In this study, we analyse collisionless simulations of clusters 
of galaxies that are generated by the technique of 'zooming 
in' on a region of interest (Tormen et al. 1997). In a first step, 
a cosmological simulation with sufficiently large volume is 
used to allow the selection of a suitable target cluster. For 
this purpose, we employed the GIF- ACDMt model carried 
out by the Virgo consortium. It has cosmological parameters 
Qo = 0.3, Q.a = 0.7, h = 0.7$, spectral shape T = 0.21, and 
was cluster-normalized to as = 0.9. The simulation followed 
256 3 particles of mass 1.4 x 10 10 /i _1 Mo within a comovmg 
box of size 141.3 /i _1 Mpc on a side. Note that this simulation 
is one of the models recently studied by KCDW. We selected 
the second most massive cluster that had formed in the GIF 
simulation for further study. This cluster has a virial mass 
8.4 x 10 14 /i _1 M , and it appears to be well relaxed at the 
present time. 

In a second step, we simulated the formation of this 
cluster again, using greatly increased mass and force res- 
olution. To this end, the particles in the final GIF-cluster 
and in its immediate surroundings were traced back to their 
Lagrangian region in the initial conditions. The correspond- 
ing part of the displacement field was then sampled using a 
glass-like particle distribution with smaller particle masses 
than in the GIF simulation. Due to the increase in reso- 
lution, the fluctuation spectrum could now be extended to 
smaller scales. We added a random realization of this ad- 
ditional small-scale power, while we kept all the waves on 
larger scales that had been used in the GIF simulation. 

Outside this central high-resolution region, we gradually 
degraded the resolution by using particles with masses that 
grow with distance from the center. In this 'boundary re- 
gion', we employed a spherical grid whose spacing grew with 
distance from the high resolution zone. The spherical bound- 
ary region extends to a total diameter of 141.3 h~ 1 Mpc, 
which is just the box size of the original GIF simulation. Be- 
yond this region, we assumed vacuum boundary conditions, 
i.e. a vanishing density fluctuation field. Using comoving co- 
ordinates, we then evolved the simulations to redshift z = 
with our new parallel tree-code GADGET. This code uses 
individual timesteps for all particles, and was designed to 
run on massively parallel supercomputers with distributed 

I The GIF project is a joint effort by astrophysicists in Germany 
and Israel. 

t We employ the convention Hq = 100?ikms _1 Mpc -1 . 
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Table 1. Numerical parameters of our cluster simulations. All four simulations compute the evolution of the same cluster, assuming a 
ACDM universe with cosmological parameters f2o = 0.3, (i\ = 0.7, V = 0.21, <T8 = 0.9, and h = 0.7. The simulations follow a sphere of 
matter with comoving diameter 141h~ 1 Mpc. In the Table, we give the particle mass m p used in the central high-resolution zone, the 
starting redshift z atar t, the gravitational softening e in the high-resolution zone, the number N hl of high-resolution particles, the number 
Abnd of boundary particles, the total number iVtot of particles, and the number N p of processors used in each of the simulations S1-S4. 
The gravitational softening was kept fixed at the given values in physical coordinates below redshift 2 = 9, and in comoving coordinates 
above this redshift. 
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Figure 1. The projected mass density fields in slices of thickness 10/i -1 Mpc around the cluster center in the original GIF simulation 
(left), and in the S3 resimulation (right). The left image is 141h _1 Mpc on a side, and the white square marks the region (85/i -1 Mpc 
on a side) that is displayed in the image of the resimulation on the right. In the right panel, you may notice small traces of the spherical 
grid used to represent the density field in the boundary region. Note that these residuals of the grid structure are just seen because of 
projection effects that arise in the visualization technique. 



memory. Parallelization is achieved explicitly using the com- 
munication library of the Message Passing Interface (MPI). 
A full account of the numerical and algorithmic details of 
GADGET is given elsewhere (Springel et al. 2000). 

To be able to study systematic effects arising from nu- 
merical resolution, we simulated the same cluster several 
times, increasing the resolution step by step. In the first 
step, refered to as simulation 'SI' from here on, the particle 
mass was just about two times smaller than in the original 
GIF simulation. For the high-resolution zone of SI, we used 
a total of 450000 particles and a gravitational softening of 
e = 6.0/i -1 kpc. The boundary region was represented with 
an additional ~ 3 million particles. This relatively high num- 
ber of boundary particles was chosen with the sequence of 



our planned simulations in mind. Except for slight changes 
at its inner rim, we kept the sampling of the boundary re- 
gion fixed for the other simulations, where we populated the 
central zone with many more particles. In simulation 'S2', 
we used 2 million high-resolution particles, and we reduced 
the gravitational softening to 3.0 /i -1 kpc. In simulation 'S3', 
we employed a total of 13 million high-resolution particles 
with mass 2.4 x 10 8 /i _1 M0, and a softening of 1.4/i _1 kpc. 
Finally, in simulation 'S4' we used 66 million particles for 
the high-resolution zone, pushing the particle mass down 
to 4.68 x 10 7 Ii~ 1 Mq and the spatial resolution down to 
0.7/i _1 kpc. In each case, roughly one third of the particles 
in the high-resolution zone ends up in the virialized region of 
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Figure 2. The dark matter distribution of the S4 cluster at z = 0. The image shows all the mass in a box of sidelength 4 ft -1 Mpc around 
the cluster center. To render the substructure visible, particles have been weighted by their local density (computed by adaptive kernel 
estimation), and a logarithmic color scale has been applied. Note that the small bright dots that are visible in the cluster should not be 
mistaken as noise - they are in fact self-bound subhaloes and correspond to surviving cores of haloes that have fallen into the cluster at 
some earlier time. 



the final cluster. This means that S4 resolves a single object 
with about 20 million particles. 

In Table 1, we summarize important numerical param- 
eters of our simulations. Note that we have softened gravity 
using a spline kernel. Our cited values for e are such that 
the gravitational potential of a point mass at zero lag is 
$ = — Gm/f, and that the softened force becomes Newto- 
nian at a distance 2.8 e. We have kept the softening length 
fixed in physical coordinates below redshift z = 9, and in 
comoving coordinates at higher redshift. The softening for 
the boundary particles was set to much larger values, in 
an inner shell around the cluster to 15 h~ kpc, and further 
outside to 75/i _1 kpc. For all four simulations, we stored 51 
outputs, logarithmically spaced in expansion factor between 
redshifts z = 20 and z = 0. 

© 0000 RAS, MNRAS 000, 000-000 



In Figure 1, we show two images comparing the original 
density field of the GIF-simulation with that of our S3 res- 
imulation. The filaments of dark matter around the cluster 
are nicely reproduced by S3, even relatively far away from 
the cluster, where the resolution of S3 has already fallen be- 
low that of the GIF simulation. It is also apparent that S3 
has much higher resolution in the cluster itself. 

However, the vast increase in resolution that this new 
set of simulations offers is perhaps best appreciated if we 
zoom in onto the cluster directly. In Figure 2, we show an 
image of the dark matter distribution in a box of side-length 
4/i _1 Mpc, centered on the cluster that formed in S4, the 
simulation with our highest resolution. Note that there is 
essentially no particle noise visible in this picture; all the 
small bright features are genuine self-gravitating subhaloes. 
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So dark matter haloes forming in CDM cosmologies exhibit 
a remarkable richness of substructure, and are thus quite far 
from the smooth, over-merged objects suggested by numer- 
ical work on cluster formation until a few years ago. This 
new view of CDM dynamics has only recently been fully es- 
tablished, with some of the most important work done by 
Ghigna et al. (1998); Moore et al. (1999) and Klypin et al. 
(1999). Our best simulation achieves an even larger dynamic 
range than previous work; as we will see, the virialized re- 
gion of the S4 cluster contains nearly 5000 self-gravitating 
subhaloes. 



3 MODELING GALAXY FORMATION USING 
N-BODY MERGING TREES 

In the following, we briefly summarize our specific imple- 
mentation of the techniques developed by KCDW to com- 
bine semi-analytic models for galaxy formation with dark 
matter merging history trees constructed directly from cos- 
mological N-body simulations. We will later extend this for- 
malism to include dark matter substructure, and we will be 
especially interested in any changes of the results arising 
from that. 

There are essentially two main parts in the modeling: 
(1) The measurement of dark matter merging trees from a 
sequence of simulation outputs. (2) The implementation of 
the actual semi-analytic recipes for galaxy formation on top 
of these merging trees. Both parts of the modeling are tech- 
nically complex and warrant a detailed discussion, which we 
provide in the following sections for the sake of complete- 
ness. Readers who are primarily interested in the results of 
our modelling may want to proceed directly to Section 5 
upon a first reading of this paper. 

In this Section, we start by describing our implemen- 
tation of the techniques of KCDW. First we treat the con- 
struction of the merging trees, then the physics of galaxy 
formation. In Section 4, we describe what we change in the 
two parts to allow the inclusion of subhalo information. 



3.1 Following the merging trees 

For each simulation output, we compile a list of dark matter 
haloes with the friends-of-friends (FOF) algorithm using a 
linking length of 0.2 in units of the mean interparticle sepa- 
ration. We only include groups with at least 10 particles in 
the halo catalogue. The majority of such haloes are already 
stable, i.e. particles found in a 10-particle group at one out- 
put time are almost all part of the same halo in subsequent 
simulation outputs. For each halo, we also determine the 
most-bound particle within the group, where 'most-bound' 
here refers to the particle with the minimum binding energy. 

We then follow the merger tree of the dark matter from 
output to output. A halo Hb at redshift zb is defined to be 
a progenitor of a halo Ha at redshift za < zb, if at least 
half of the particles of Hb are contained within Ha, and the 
most bound particle of Hb is contained in Ha, too. These 
definitions already suffice to uniquely define the dark matter 
merging trees. 



3.1.1 Defining a galaxy population 

So far we are just dealing with catalogues of dark mat- 
ter haloes. We now supplement this with the notion of a 
galaxy population with physical properties given by the semi- 
analytic techniques. In our formalism, each dark halo carries 
exactly one central galaxy, and its position is given by the 
most-bound particle of the halo. Only the central galaxy is 
supplied with additional gas that cools within the halo. 

A halo can also have one or several satellite galaxies, 
where the position of each of them is given by one of the 
particles of the halo. Satellites are galaxies that had been 
central galaxies themselves in the past, but their haloes have 
merged at some previous time with the larger halo they now 
reside in. Satellite galaxies orbit in their halo and are as- 
sumed to merge with the central galaxy on a dynamical 
friction timescale. Note that they are cut off from the sup- 
ply of fresh cool gas, so they may only form stars until their 
internal reservoir of cold gas is exhausted. 

Finally, we define a class of field galaxies, which arc 
introduced to keep track of satellites whose particles are 
currently not attached to any halo, for example because they 
have been ejected out of their parent halo. Usually, these 
"lost" field galaxies are accreted onto a halo later on. 

3.1.2 Following the galaxy population in time 

At a given output time, we therefore deal with a galaxy pop- 
ulation consisting of central galaxies, satellite galaxies, and 
field galaxies, each attached to the position of a simulation 
particle. Starting at the first output time at high redshift 
(when the first haloes have formed) , we initialize the galaxy 
population with a set of central galaxies, one for each halo, 
with stellar mass, cold gas mass, and luminosity set to zero. 
The physical properties of these galaxies are then evolved to 
the next output time, where we obtain a new galaxy popu- 
lation based on a combination of semi-analytic prescriptions 
and the merging history of the dark matter. Propagating 
this scheme forward in time from output to output we ob- 
tain the galaxy population at the present time, and at all 
output times at higher redshift. 

We now describe in more detail our rules and prescrip- 
tions for this evolution. Beginning with the galaxy popula- 
tion at redshift zb , we first generate the galaxies of the new 
population at redshift za < zb based on the merging history 
of the dark matter. Using the group catalogues of the corre- 
sponding simulation outputs and the galaxy population at 
zb, we construct an 'initial' set of galaxies at za as follows: 

(i) Each galaxy at zb is assigned to its new halo at za- 
If the particle used to tag a galaxy does not reside in any 
halo, the galaxy becomes a field galaxy. 

(ii) Each halo at za selects its central galaxy as the cen- 
tral galaxy of its most massive progenitor. This central 
galaxy is repositioned to the position of the most-bound par- 
ticle, i.e. a new particle tagging the galaxy is selected. The 
central galaxies of all other progenitors become satellites of 
the halo. 

(iii) If a halo has no progenitors, a new central galaxy is 
created at the position of its most-bound particle. In the 
event that the halo contains one or more galaxies (particles 
recovered from the field) , the central galaxy is picked as the 
most massive of these. 
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Once the new set of galaxies is generated in this way, 
the properties of the galaxies are evolved according to the 
physical prescriptions described below, resulting finally in 
the new galaxy population at redshift za- Note that some of 
the satellite galaxies generated in the initial set for za will 
merge with central galaxies during this evolution, and thus 
not necessarily be part of the 'final' population at za- 



3.2 Physical evolution of the galaxy population 

We model the following physical processes: (1) Radiative 
cooling of hot gas onto central galaxies. (2) Transformation 
of cold gas into luminous stars by star formation. (3) Reheat- 
ing of cold gas, or its ejection out of the halo, by supernova 
feedback. (4) Orbital decay of satellites and their merging 
with central galaxies. (5) Spectrophotometric evolution of 
the luminous stars. (6) Simplified morphological evolution 
of galaxies. Below we detail the physical parameterizations 
adopted for these processes. 



3.2.1 Gas cooling 

Gas cooling is modeled as in White & Frenk (1991). We 
assume that the hot gas within a dark halo is distributed 
like an isothermal sphere with density profile p g (r). Then 
the local cooling time t coo \(r) can be defined as the ratio 
of the specific thermal energy content of the gas, and the 
cooling rate per unit volume, viz. 



fcTp g (r) 



2 Jlm v nl{r) A(T, Z) ' 



(1) 



Here JLm p is the mean particle mass, n e (r) the electron den- 
sity, and A(T, Z) the cooling rate. The latter depends quite 
strongly on the metallicity Z of the gas, and on the virial 
temperature T = 35.9 (Vvir/kms -1 ) 2 K of the halo. We em- 
ploy the cooling functions computed by Sutherland & Do- 
pita (1993) for collisional ionisation equilibrium to represent 
A(T, Z), but we restrict ourselves to primordial metallicity 
in the present study. 

We define the cooling radius 

T cool tlS the radius for which 
t CO oi is equal to the time for which the halo has been able to 
cool 'quasi-statically'. We will approximate this time with 
the dynamical time R v i Y /V v i Y of the halo. If the cooling ra- 
dius lies well within the virial radius R V [ T of a given halo, we 
then take the cooling rate to be 



dM c , 



cool 



(2) 



dt dt 

Note that we define the virial radius 7? v ir of a FOF-halo 
as the radius of a sphere which is centered on the most- 
bound particle of the group and has an overdensity 200 with 
respect to the critical density. We take the enclosed mass 
Mvir = 100-ff 2 i?vir/G as the virial mass, and we define the 
virial velocity as V} 1T = GM vlI / R vlI . 

Adopting an isothermal sphere for the distribution of 
the hot gas of mass Mhot within the halo, i.e. 

M hot 



the cooling rate is then given by 
dM cc 



M ho t r coo i 



dt 



-Rvir 2£ CO ol 



(3) 



(4) 



At early times or for low-mass haloes the cooling radius can 
be much larger than the virial radius. In this case, the hot 
gas is never expected to be in hydrostatic equilibrium, and 
the cooling rate will essentially be limited by the accretion 
rate onto the central galaxy. We approximate this rate with 



dM a , 



dt 



M hot 

2^cool 



(5) 



which corresponds to r coo i = R viT in equation (4). We adopt 
the minimum of equations (4) and (5) as our actual cool- 
ing rate. As in KCDW we find it necessary to implement 
an upper limit on the circular velocity of haloes in which 
cooling gas is allowed to settle onto a central galaxy. This 
accounts for the fact that observed cluster cooling flows do 
not form stars at the observed apparent cooling rate (which 
corresponds to that estimated above). We follow KCDW in 
setting this upper limit at a circular velocity of 350kms _1 . 



3. 2. 2 Star formation 

In this work, we model the star formation rate of a galaxy 
as 



dM* M C oid 

— ; — = a 

dt 



tdy 



(6) 



where M co id is the mass of its cold gas, and tdyn is the dy- 
namical time of the galaxy. We approximate the latter as 



tdy 



geff 
Kir' 



(7) 



with i? e ff = 0.17? v ir, i.e. we set the effective stellar radius to 
be a fixed fraction of the virial radius. Notice that at a fixed 
redshift, R v i T is proportional to VW, hence td yn depends only 
on redshift. The dimensionless parameter a regulates the ef- 
ficiency of star formation and is treated as a free parameter. 
In this work, we keep a constant in time, but we remark 
that a redshift dependence of a may be required to provide 
a better understanding of the rapid evolution of the number 
density of luminous quasars (Kauffmann & Haehnelt 2000) 
and to match the observed abundance of Lyman break galax- 
ies at z ~ 3 (Somerville et al. 2000). Once a galaxy falls into 
a larger halo and becomes a satellite, the values of i? v ; r and 
Vvir are not changed any more. The galaxy can then continue 
to form stars until its reservoir of cold gas is exhausted, but 
it does not receive new cold gas by cooling processes. 



3.2.3 Feedback 

Assuming a universal initial mass function (IMF) , the energy 
released by supernovae per formed solar mass is rysN-EsN, 
where ??sn gives the expected number of supernovae per 
formed stellar mass, and Esn is the energy released by each 
supernova. The formation of a group of stars with mass AM* 
will thus be accompanied by the release of a feedback en- 
ergy of ^sn-EsnAM*, where we adopt tisk = 5.0 x 10~ 3 Mq 1 , 
based on the Scalo (1986) IMF, and E S n = 10 51 erg. 

One major uncertainty is how this energy affects the 
evolution of the interstellar medium, and how the star for- 
mation rate is regulated by it (Springel 2000). We here as- 
sume that the feedback energy reheats some of the cold gas 
back to the virial temperature of the dark halo. The amount 
of gas reheated by this process is then 
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AMrohcat = - e — -~2 — AM*, 
3 y^r- 



(8) 



where the dimensionless parameter e describes the efficiency 
of this process. 

Following KCDW, we consider two alternative schemes 
for the fate of the reheated gas. In the retention scheme, 
the reheated gas is simply transferred from the cold phase 
back to the hot gaseous halo, and the reheated gas thus 
stays within the halo. Alternatively, in the ejection scheme 
we assume that the gas leaves the halo, and it is only re- 
incorporated into the halo at some later time. If AM c j cc 
is the total gas mass ejected by a galaxy, we model this 
reincorporation by decreasing AM c j oc to zero again on the 
dynamical timescale of the halo. 



3.2.4 Mergers of galaxies 

In CDM universes, large haloes form by mergers of smaller 
haloes. As a consequence, mergers of galaxies are an in- 
evitable process. We assume that the satellite galaxies orbit- 
ing within a dark matter halo experience dynamical friction 
and will eventually merge with the central galaxy of the 
halo. In principle, mergers between two satellite galaxies are 
also possible. These events are expected to be rare, but they 
do happen occasionally as we show in the companion paper. 
For the moment, we neglect these events but we will take 
them into account in our subhalo-scheme later on. 

N-body simulations by Navarro et al. (1995) suggest 
that the merging timescale can be reasonably well approxi- 
mated by the dynamical friction timescale 



2 c 



GM sat In A 



(9) 



The formula is valid for a small satellite of mass M sat orbit- 
ing at a radius r c in an isothermal halo of circular velocity 
V c . The function /(e) describes the dependence of the de- 
cay on the eccentricity of the satellites' orbit, expressed in 
terms of e = J/J C (E), where J C (E) is the angular momen- 
tum of a circular orbit with the same energy as the satellite. 
The function /(e) is well approximated by /(e) ~ e 0,78 , for 
e > 0.02 (Lacey & Cole 1993). C is a constant with value 
C ~ 0.43, and In A is the Coulomb logarithm. 

We follow KCDW and approximate r c with the virial 
radius of the halo when the satellite first falls into it. To de- 
scribe the orbital distribution, we adopt the average value 
(/(e)) ~ 0.5, computed by Tormen (1997). Note that this 
differs slightly from KCDW who drew a random orbit uni- 
formly from e € [0.02, 1]. We identify the mass of the satel- 
lite with the virial mass of the galaxy at the time when it 
was last a central galaxy, and we approximate the Coulomb 
logarithm with In A = (1 + M viT /M BS ,t). 

When a small satellite merges with a central galaxy, we 
transfer all its stellar mass to the bulge component of the 
central galaxy, and we update the photometric properties of 
this galaxy accordingly. Similarly, the cold gas of the satel- 
lite is transferred to the disk of the central galaxy. If the 
mass ratio between the stellar components of the merging 
galaxies is larger than some threshold value (we adopt 0.3 
for that), the merger destroys the disk of the central galaxy 
completely, and all stars form a single spheroid, i.e. they 
generate a bulge. This is called a major merger. In addition, 



we assume that all the cold gas left in the two merging galax- 
ies is rapidly consumed in a starburst. The stars created in 
this burst are also added to the bulge component. Since the 
central galaxy is fed by a cooling flow, it can grow a new 
disk component later on. 



3.2.5 Spectrophotometric evolution 

Photometric properties of our model galaxies can be con- 
structed using stellar population synthesis models (Bruzual 
& Chariot 1993). In these models, the number of stars that 
initially form in each mass range is computed according to 
the initial mass function. The stars then evolve along the- 
oretical evolutionary tracks. In this way, the spectra and 
colors of a stellar population formed in a short burst of star 
formation can be followed as a function of time. Once the 
evolution F v (t) of the spectral energy distribution (SED) of 
a single age population of stars is known, the SED S v (t) of 
a galaxy can be computed as 



S v (t) = [ F v (t-t')M,(t')dt' 
Jo 



(10) 



from its star formation history M*(i). Upon convolution 
with standard filters, colors and luminosities in the desired 
bands can be obtained. In principle, this technique also al- 
lows the redshifting of spectra, and the incorporation of k- 
corrections to make direct contact with observational pho- 
tometric data at high redshift. In this work we use updated 
evolutionary synthesis models by Bruzual & Chariot (in 
preparation), which have been computed for solar metallic- 
ity. Note that in this paper we do not attempt to model the 
effects of dust. These effects can be be quite substantial in 
late-type galaxies and would significantly affect a number of 
our results (c.f. Somerville & Primack 1999). In particular, 
when normalising our models to the observed Tully-Fisher 
relation, the inclusion of dust would cause us to assign a 
higher stellar mass to a galaxy of given circular velocity 

3.2.6 Morphological evolution 

Simien & de Vaucouleurs (1986) find a good correlation 
between the B-band bulge-to-disk ratio, and the Hubble- 
type T of galaxies. For a magnitude difference AM = 
Afbuigc — Mtotai they find a mean relation 

(AM) = 0.324(T + 5) - 0.054(T + 5) 2 + 0.0047(T + 5) 3 .(11) 

Following previous semi-analytic studies, we assign mor- 
phologies based on this equation. Specifically, we will usu- 
ally classify galaxies with T < —2.5 as ellipticals, those with 
-2.5 < T < 0.92 as SO's, and those with T > 0.92 as spirals 
and irregulars. However, we may allow a shift in the bound- 
aries between the three classes to obtain a better match of 
our models with the observed relative abundances of these 
three morphological types (see section 5.3). Note that galax- 
ies without any bulge are classified as type T = 9. 



3.3 More implementation details 

Our practical implementation of the physical evolution of 
the galaxy population includes the following steps. We first 
estimate merger timescales for those satellites that have 
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newly entered a given halo, i.e. the galaxies that had not 
been contained in the largest progenitor of the halo. This 
'merger clock' is then decreased with time in the subsequent 
evolution, and the satellite will be merged with the central 
galaxy when this time has elapsed. Note that the merger 
clock may be reset before the merger happens if the halo 
containing the satellite merges with a larger system. 

We then compute the total amount of hot gas available 
for cooling in each halo. Assuming on average a universal 
abundance of baryons equal to the primordial one, this is 
simply 



M hot = f h M v 



Mf> + M?> + M 



(i) 



r(i) 



(12) 



where the sum extends over all galaxies within the halo. Here 
/ b = fib /Ho denotes the baryon fraction of the universe. Us- 
ing the cooling model of Section 3.2.1, we then estimate the 
cooling rate onto each central galaxy, and we keep this rate 
constant during the time AT between the two simulation 
outputs. 

Once these quantities are known, we solve the sim- 
ple differential equations describing star formation, cooling 
and feedback. We typically use a number of N ~ 50 small 
timesteps of size At — AT /N for this purpose. At each 
of these small steps, new cold gas is added to the central 
galaxies. For each galaxy, we form some stellar mass AM* 
according to its star formation rate, and we update its cur- 
rent and future photometric properties accordingly. The cold 
gas mass of each galaxy is reduced by the amount of stars 
formed, and by the mass of the gas that is reheated or ejected 
by supernova feedback. 

At the end of each of the small steps, the merger clocks 
of the satellites are reduced by At. If a satellite's merging 
time falls below zero, it is merged with the central galaxy of 
its parent halo. In practice, this means that the luminosity, 
the stellar mass, and the gas mass of the satellite are trans- 
fered to the central galaxy, and that the satellite is removed 
from the list of galaxies. In addition, in the event of a major 
merger all cold gas of the central galaxy is consumed in a 
short starburst, and all stellar material is transformed into 
a spheroid. 

3.4 Choice of model parameters 

Following KCDW, we use the /-band Tully-Fisher relation to 
normalize our models, i.e. to set the free parameters a and 
e which specify the efficiency of star formation and feed- 
back, respectively. For that purpose, we consider the cen- 
tral galaxies of haloes in the periphery of the cluster, with 
morphological types corresponding to Sb/Sc galaxies. Note 
that we only use haloes that are not contaminated by heav- 
ier boundary particles. The remaining number of galaxies is 
sufficiently large to construct a well defined Tully-Fisher di- 
agram. We try to fit the velocity based /-band Tully-Fisher 
relation 



Mi - 5 log h = -21.00 - 7.68(log W - 2.5) 



(13) 



measured by Giovanelli et al. (1997). We set the velocity 
width W as twice the circular velocity, and we assume that 
the circular velocity V c of a spiral galaxy is ~ 15% larger 
than the virial velocity of that galaxies halo. This is moti- 
vated by detailed models for the structure of disk galaxies 



by Mo et al. (1998) embedded in cold dark matter haloes 
with the universal NFW profile (Navarro et al. 1996, 1997). 

Keeping other parameters fixed, we find that varying e 
changes both the slope and the zero-point of the Tully-Fisher 
relation strongly. In particular, making feedback stronger 
tilts the Tully-Fisher relation towards steeper slopes. On the 
other hand, the star formation efficiency a only very weakly 
affects the zero-point, but it has a strong effect on the gas 
mass fraction left in galaxies at the present time. 

In principle, it should be possible to specify the param- 
eters a and e using the slope and zero-point of the Tully- 
Fisher relation alone. However, the weak dependence of the 
Tully-Fisher relation on a makes this impractical. As in 
KCDW, we instead use an additional criterion and require 
that the cold gas (HI plus molecular) mass in a 'Milky- 
Way' galaxy of circular velocity V c = 220kms~ is about 
8 x 1O 9 /i _1 M . 

Note that the baryon fraction /b can strongly influence 
the cooling rates, and thus the absolute normalization of 
the models. As White et al. (1993) have shown, the baryon 
content of rich clusters of galaxies argues for a baryon frac- 
tion as high as /b = 0.1 — 0.2. This is inconsistent with big 
bang nucleosynthesis (BBN) constraints in a critical density 
universe, but can be accommodated within low-density cos- 
mologies, like the one considered in our cluster models. We 
will assume / b = 0.15 in this study, which is consistent with 
current BBN constraints. Our resulting parameter values are 
listed in Table 2. We expect that slightly different values of 
/b will produce similar results when the parameters a and 
j3 are adjusted to compensate. 



4 FOLLOWING HALO SUBSTRUCTURE 

4.1 Identification of substructure 

A basic step in the analysis of cosmological simulations is the 
identification of virialized particle groups, which specify the 
sites where luminous galaxies form. Perhaps the most popu- 
lar technique employed for this task is the friends-of-friends 
(FOF) algorithm. It places any two particles with a sepa- 
ration less than some linking length b into the same group. 
In this way, particle groups are formed that correspond to 
regions approximately enclosed by isodensity contours with 
threshold value p oc 1/b 3 . For an appropriate choice of b, 
groups are selected that are close to the virial overdensity 
predicted by the spherical collapse model. FOF is both sim- 
ple and efficient, and its group catalogues agree quite well 
with the predictions of Press-Schechter theory. 

However, FOF has a tendency to link independent 
structures across feeble particle bridges occasionally, and in 
its standard form with a linking length of b ~ 0.2 it is not 
capable of detecting substructure inside larger virialized ob- 
jects. Using sufficiently high mass resolution, recent stud- 
ies (Tormen 1997; Tormen et al. 1998; Ghigna et al. 1998; 
Klypin et al. 1999; Moore et al. 1999) were able to demon- 
strate that substructure in dense environments like groups 
or clusters may survive for a long time. The cores of the dark 
haloes of galaxies that fall into a cluster will thus remain in- 
tact, and orbit as self-gravitating objects in the smooth dark 
matter background of the cluster. In previous simulations, 
haloes falling into clusters usually evaporated quickly, and 
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the clusters exhibited little signs of substructure (e.g. Frenk 
et al. 1996). It now appears, that sufficient numerical force 
and mass resolution is enough to resolve this "overmerging" 
problem. 

The identification of substructure within dark matter 
haloes is a challenging technical problem, and several al- 
gorithms to find "haloes within haloes" have been pro- 
posed. In hierarchical friends-of-friends (HFOF) algorithms 
(Gottlober et al. 1998; Klypin et al. 1999) the linking length 
of plain FOF is reduced in a sequence of discrete steps, thus 
selecting groups of higher and higher overdensity and even- 
tually capturing true substructure. 

Clearly, the need for a well-posed physical definition of 
"substructure" arises early on in such an analysis. Most au- 
thors have required subhaloes to be locally overdense and 
self-bound. We will also adopt this requirement. Note that 
this implies that any locally overdense region within a dense 
background needs to be treated with an unbinding proce- 
dure. This is because a small halo within a larger system 
represents only a relatively small fluctuation in density, and 
a substantial amount of mass within the overdense region 
will just stream through and not be gravitationally bound 
to the substructure itself. 

Groupfinding techniques that use some criterion of self- 
boundedness include the bound density maximum (BDM) 
algorithm (Klypin et al. 1999), where the bound subset of 
particles is evaluated iteratively in spheres around a local 
density maximum. In the method of Tormen et al. (1998), 
previous simulation outputs are used to track the infall of 
particle groups into larger systems. Once such a particle 
group from the field was accreted by a cluster, they simply 
determined the subset of those particles that still remained 
self- bound. 

Another approach is followed in DENMAX (Gelb & 
Bertschinger 1994) and its offspring SKID, where particles 
are moved along the local gradient in density towards a local 
density maximum. Particles ending up in the 'same' maxi- 
mum are then linked together as a group using FOF. SKID 
has been employed by Ghigna et al. (1998) to find substruc- 
ture in a rich cluster of galaxies, and to study the statistical 
properties of the detected subgroups. 

Integrating the gradient of the density field and moving 
the particles is not without technical subtleties. For example, 
a suitable stopping condition is needed. The new algorithm 
HOP of Eisenstein & Hut (1998) tries to avoid these diffi- 
culties by restricting the group search to the set of original 
particle positions, just like FOF does. In HOP, one first ob- 
tains an estimate of the local density for each particle, and 
then attaches it to its densest neighbour. In this way a set 
of disjoint particle groups are formed. However, a number of 
additional rules are needed to link and prune some of these 
groups. For example, HOP may split up a single virialized 
clump into several pieces of unphysical shape, which have to 
be joined using auxiliary criteria. 

It appears that all of these techniques have different 
strengths and weaknesses, and that none is completely sat- 
isfactory (for example, DENMAX and HOP do not require 
the identified substructure to be gravitationally self-bound). 
We have therefore come up with a new algorithm to detect 
substructure in dark matter haloes that incorprates ideas 
from SKID, HOP, FOF, and IsoDen (Pfitzner et al. 1997), 



as well as adding some new ones. For easier reference, we 
dub this algorithm SUBFIND (for subhalo finder). 

4.2 The algorithm SUBFIND 

Our objective with SUBFIND is to be able to extract sub- 
structure, which we define as locally overdense, self-bound 
particle groups within a larger parent group. The parent 
group will be a particle group pre-selected with a standard 
FOF linking length, although SUBFIND could operate on ar- 
bitrary particle groups, or with slight modifications on all of 
the particles in a simulation at once. The use of FOF-groups 
as input data provides a convenient means to organize the 
groups according to a simple two stage hierarchy consisting 
of 'background group' and 'substructure'. 

Note that it is unlikely that we lose any substructure 
by restricting the search to ordinary FOF-groups. FOF may 
accidentally link two structures, a case SUBFIND will be 
able to deal with, but we rarely expect FOF (with a linking 
length of 0.2) to split a physical structure into two parts. 

In SUBFIND, we begin by computing a local estimate 
of the density at the positions of all particles in the input 
group. This is done in the usual SPH-fashion, i.e. the local 
smoothing scale is set to the distance of the iVdcns nearest 
neighbour, and the density is estimated by kernel interpola- 
tion over these neighbours. The particles may be viewed as 
tracers of the three-dimensional dark matter density field. 
We consider any locally overdense region within this field 
to be a substructure candidate. More specifically, we define 
such a region as being enclosed by an isodensity contour that 
traverses a saddle point. How can one find these regions? 
Imagine lowering a global density threshold slowly within 
the density field. For most of the time, isolated overdense 
regions will just grow in size as the threshold is lowered, ex- 
cept for the moments when two separate regions coalesce to 
form a common domain. Note that at these instances, the 
contours of two separate regions join at a saddle point. As 
a result the topology of the isodensity contour changes as 
well. 

Our algorithm tries to identify all locally overdense re- 
gions by imitating such a lowering of a global density thresh- 
old. To this end, we sort the finite number of particles ac- 
cording to their density, and we 'rebuild' the particle distri- 
bution by adding them in the order of decreasing density. 
Whenever a new particle i with density pi is considered, we 
find the N ngh nearest neighbors within the full particle set. 
Within this set Ai of N ng b particles, we also determine the 
subset of particles with density larger than pi, and among 
them we select a set Bi containing the two closest particles. 
Note that this set may contain only one particle, or it may 
be empty. We now consider three cases: 

(i) The set Bi is empty, i.e. among the iV ngb neighbors 
is no particle that has a higher density than particle i. In 
this case, particle i is considered to mark a local density 
maximum, and it starts growing a new subgroup around it. 

(ii) If Bi contains a single particle, or two particles that 
are attached to the same subgroup, the particle i is also 
attached to this subgroup. 

(iii) Bi contains two particles that are currently attached 
to different subgroups. In this case, the particle i is consid- 
ered to be a saddle point, and the two subgroups labeled 
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Figure 3. Example for a subhalo identification with SUBFIND. The top left panel shows a small FOF-group (44800 particles), identified 
at z = in the vicinity of the S2 cluster. SUBFIND identifies 56 subhaloes within this group, the largest one forms the background halo 
and is shown on the top right, while the other 55 subhaloes are plotted on a common panel on the lower left. In this example, the total 
mass in all the "true" subhaloes 2-56 is about 8% of the group mass. Particles not bound to any of the subhaloes form "fuzz" , and are 
displayed on the lower right. These particles primarily lie close to the outer edge of the group. Spatial coordinates are given in h~ 1 kpc. 



by the particles in Bi are registered as subhalo candidates. 
Afterwards, the particle i is added by joining the two sub- 
groups to form a single subgroup. Note that all subhalo can- 
didates will be examined for self-boundedness later on in the 
algorithm. 

Working through this scheme results in a list of subhalo 
candidates, which can be efficiently stored in a suitably ar- 
ranged link-list structure. Note that a given particle can be 
member of several different subhalo candidates, and that the 



algorithm is in principle fully capable of detecting arbitrary 
levels of "subhaloes within subhaloes" . 

Up to this point, the construction of subhalo candidates 
has been based on the spatial distribution of particles alone. 
A more physical definition of substructure is obtained by 
adding the requirement of self-boundedness. We therefore 
subject each subhalo candidate to an unbinding procedure 
to obtain the "true" substructure. To this end, we succes- 
sively eliminate particles with positive total energy, until 
only bound particles remain. We perform the unbinding in 
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physical coordinates, where we define the subhalo's center 
as the position of the most bound particle, and the velocity 
center as the mean velocity of the particles in the group. We 
then obtain physical velocities with respect to this group 
center by adding the Hubble flow to the peculiar velocities. 
Finally, if more than a minimum number of iVngb particles 
survive the unbinding, we refer to these particles as a sub- 
halo. 

An important issue remains of how one should deal with 
complications arising from the assignment of particles to sev- 
eral different subhalo candidates. This does not only occur if 
one deals with genuine "substructure within substructure" , 
but is actually quite typical for the algorithm. For exam- 
ple, imagine a large halo containing several small subhaloes. 
Whenever one of the small haloes 'separates' from the main 
halo, two subhalo candidates are generated according to case 
(iii) of the algorithm. Each time the larger of these groups 
describes the bulk of the main halo, which will thus appear 
several times as a subhalo candidate, although one would like 
to consider it only once as an independent physical struc- 
ture. 

We approach this issue by considering only the smaller 
subhalo candidate at each branch of the tree generated by 
the saddle points. This is based on the notion that we want 
to examine substructure within some larger object, and this 
'background' object is expected to have larger mass than 
the actual substructure. In addition, we process the subhalo 
candidates in the inverse sequence as they have been gen- 
erated, i.e. we work through the saddle points from low to 
high density. In this way, a smaller subhalo within a larger 
subhalo will always be processed later than its parent sub- 
halo. As we consider the subhalo candidates in this order, 
each particle carries a label indicating the subhalo it was 
last detected to reside in. If in the process it is found to be 
contained also within a smaller subhalo, this label will be 
overwritten by the new subgroup identifier. 

In this way, the complexity of our analysis is reduced by 
assigning each particle at most to one subhalo. We are still 
able to detect a hierarchy of small subhaloes within larger 
subhaloes, albeit at the expense of reducing the latter by the 
particles contained at deeper levels of the hierarchy. How- 
ever, this usually does not affect the corresponding parent 
subhalo strongly, since the mass of any substructure within 
a larger group is usually small compared to that of the par- 
ent group. Nevertheless, it can happen that the extraction of 
subhaloes unbinds some additional particles from the parent 
subhalo. For this reason, we check all the disjoint subhaloes 
at the end of the process yet again for self-boundedness. 
Here, we also assign all particles not yet bound to any sub- 
halo to the "background halo" of the group, which we define 
as the largest subhalo within the original FOF input group, 
and we check whether they are at least bound to this struc- 
ture. If not, these particles represent 'fuzz', identified by 
FOF to belong to the group, but not (yet) gravitationally 
bound to it. 

In summary, SUBFIND decomposes a given particle 
group into a set of disjoint self-bound subhaloes, each iden- 
tified as a locally overdense region within the density field of 
the original structure. The algorithm is spatially fully adap- 
tive, and it has only two free parameters, iV dcns and N ngb . 
The latter of these parameters sets the desired mass res- 
olution of structure identification, and we usually employ 



N ns b = 10 for this purpose. The results are quite insensitive 
to the other parameter, the number iVdens of SPH smoothing 
neighbours, which we typically set to a value slightly larger 
than iV n gb. 

Finally, we note that any efficient practical implemen- 
tation of the algorithm requires the use of hierarchical tree- 
data structures, and fast techniques to find nearest neigh- 
bours and gravitational potentials. For this purpose, we em- 
ploy techniques borrowed from our tree-SPH code GADGET. 

In Figure 3, we show a typical example of substructure 
identified using SUBFIND. We selected a small group from 
the periphery of the S2-cluster, at z = 0, for this illustration. 
By eye, one can clearly spot substructure embedded in the 
FOF-group. The algorithm SUBFIND finds 56 subhaloes in 
this case. The largest one is the 'background' halo, shown 
in the top right panel of Fig. 3. It represents the backbone 
of the group, with all its small substructure removed. This 
substructure is made up of 55 subhaloes, which are plotted in 
a common panel on the lower left. While this example shows 
an isolated and well relaxed halo, we note that SUBFIND also 
performs well in cases where FOF links structures across 
feeble particle bridges, or when haloes are in the process of 
merging. In these cases, the algorithm reliably decomposes 
the FOF group into its constituent parts. 

4.3 Subhaloes in the SI, S2, S3, and S4 clusters 

Substructure within dark matter haloes is in itself a highly 
interesting subject that merits detailed investigation. What 
are the structural properties of subhaloes within haloes? 
What is their distribution of sizes and masses? What is their 
ultimate fate? Answers to these questions are highly relevant 
for a number of diverse topics such as galaxy formation, the 
stability of cold stellar disks embedded in dark haloes, or 
the weak lensing of galaxies in clusters (Geiger & Schneider 
1998, 1999). 

Tormen et al. (1998), Ghigna et al. (1998) and Klypin 
et al. (1999) have addressed some of these questions, and it 
will be interesting to supplement their work with results ob- 
tained from our new group-finding technique. However, such 
a study is beyond the scope of the present paper, where we 
focus on semi-analytic models for the galaxy population of 
the cluster. We therefore defer a detailed statistical analy- 
sis of substructure to future work, and just report the most 
basic properties of the subhaloes identified in the clusters at 
redshift z = 0. 

Especially interesting is a comparison of the mass spec- 
tra of subhaloes detected in the four cluster simulations SI, 
S2, S3 and S4. In this sequence of simulations, the mass and 
force resolution increases substantially, which clearly leads 
to a larger number of resolved subhaloes: 118 subhaloes are 
detect in the final cluster of SI, 496 in S2, 1848 in S3, and 
4667 in S4. 

The increase in numerical resolution thus unveils an 
enormous richness of structure in the cluster. As an exam- 
ple, we show a plot of the substructure in the S2-cluster in 
Figure 4. Note that almost all of the additional substructure 
that becomes visible with still higher numerical resolution is 
in objects of smaller mass than the ones resolved in Figure 4. 
This is seen in a comparison of the cumulative and differen- 
tial mass spectra, as shown in Figure 5. SI, S2 and S3 are 
actually well capable of resolving all the objects above their 
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Figure 4. Substructure in the S2 cluster at z = 0. The top left panel shows a color-coded projection of the FOF-group that contains 
the cluster. To highlight the substructure, particles have been given a weight proportional to the local dark matter density. In the top 
right panel we show the largest subhalo identified by SUBFIND, i.e. the background halo. The lower right shows the 495 other subhaloes 
identified in the object. Finally, on the lower right, we plot circles at the positions of each identified subhalo, with radius proportional to 
the third root of the particle number in the subhalo. Note that we actually found subhaloes within subhaloes in this example. 



respective resolution limits. Even close to their resolution 
limits they predict the right abundance of subhaloes of a 
given mass. 

4.4 Merging trees using substructure 

We now describe our methods to construct merging trees 
that take the identified subhaloes into account. Note that 
SUBFIND classifies all particles of a given FOF-group ei- 
ther as lying in a bound subhalo, or as unbound. Any or- 



dinary FOF-group lacking substructure will also appear in 
the list of subhaloes identified by SUBFIND, albeit only 
with its bound subset of particles. Since the requirement 
of self-boundedness is a reasonable physical condition for 
the definition of dark haloes, we therefore consider the list 
of haloes generated by SUBFIND as the 'source list' for our 
further analysis. Notice that we use the FOF-groups as con- 
venient 'containers' to establish a simple two stage hierarchy 
of haloes. We define the largest subhalo in a given FOF- 
group as the main halo hosting the central galaxy. All other 
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Figure 5. Subhalo mass functions in the four clusters SI, S2, S3 
and S4 at rcdshift z = 0. In the top panel, we plot the cumulative 
number N(m) of subhaloes with masses larger than m. The short 
vertical lines mark the ends of the graphs for the simulations SI 
(lowest resolution), S2, and S3 (second highest resolution). The 
agreement between the four simulations is quite good. This good 
agreement is also present in the differential mass function dN/dm, 
which we show in the lower panel. 



subhaloes found within the group will be considered to be 
substructure of this main halo. 

In the following discussion, the term subhalo will refer 
to any self-bound structure identified by SUBFIND, even if 
it is just the self-bound part of an ordinary FOF-group that 
has no real substructure. In addition, we adopt the following 
definitions: 

A subhalo Sb at redshift zb is a progenitor of a sub- 
halo Sa at redshift za < zb, if more than half of the JV link 
most-bound particles of Sb end up in Sa- This definition 
concentrates on the most-bound core of each structure, and 
we have found it to be very robust in tracing subhaloes be- 
tween different output times. One can choose AT link as some 
fraction of the number of particles of subhalo Sb, say all 
or half of them. However, such a condition may fail if Sb is 
deprived of its outer halo between two outputs, as it may oc- 
casionally happen when a structure of relatively large mass 
falls into a cluster. We have found that setting AT link = 10, 
equal to our lower particle limit for group identification, can 
satisfactorily treat even these cases. Note that our notion 
of 'most-bound' refers to the most negative binding energy. 
SUBFIND automatically stores each subhalo in the order of 
increasing binding energy to facilitate this kind of linking, 



i.e. the subhaloes are effectively stored from inside to out- 
side. 

We call a subhalo Sb at redshift zb a progenitor of a 
FOF-group Ga at redshift za < zb , if more than half of the 
particles of Sb are present in Ga • We also call a FOF-group 
Gb a progenitor of a subhalo Sa, if more than half of the 
particles of Sa are contained in Gb- Note that in this latter 
condition we deliberately used the particles of Sa to define 
'membership' in one of the FOF-group at higher redshift. 

We expect that once a self-bound structure has formed, 
most of its mass will remain in bound structures in the fu- 
ture. However, occasionally it may happen that a group that 
was just barely above our specified minimum particle num- 
ber at one output time falls below this limit at the next 
output time, for example because the group is evaporated 
by interactions with other material, or because of noise in 
the identification of groups with size close to our identifica- 
tion threshold. We call such groups volatile, and drop them 
from our analysis. More precisely: All FOF-groups without 
any bound subhalo are considered to be volatile and disre- 
garded. In addition, if a subhalo is not progenitor to any 
other subhalo, and not progenitor to any non- volatile FOF- 
group, it is considered to be volatile too, and dropped from 
further analysis. Note that basically all subhaloes eliminated 
in this way have particle numbers very close to the detection 
threshold. 

After the elimination of volatile subhaloes, we link sub- 
haloes between pairs of successive simulation outputs. By 
construction, every subhalo may have several progenitors, 
but itself can only be progenitor for at most one other sub- 
halo. In fact, due to the elimination of volatile subhaloes, a 
subhalo Sb will always be progenitor to another subhalo, or 
at least to a FOF-group Ga- If only the latter is the case, 
we treat the subhalo Sb as a progenitor of the main halo 
within the FOF-group. 

There remains then the important case that a subhalo 
Sa has no progenitor. If it also has no progenitor FOF-group, 
we call this subhalo a new galaxy, which is considered to have 
newly formed between the two output times. In the galaxy 
formation scheme, new galactic 'seeds' will be inserted into 
these subhaloes. 

If however the subhalo Sa has no progenitor subhalo, 
but a progenitor FOF-group Gb, it is likely that either the 
subhalo represents a chance fluctuation, or that it was over- 
looked in the identification process at the time zb ■ In the lat- 
ter case the corresponding structure will have been merged 
with the main halo, so we drop these subhaloes. Their abso- 
lute number and the total mass contained in them is always 
very small. 

In summary, the above identification scheme allows a 
detailed tracing of the dark matter merging history tree, 
from the past to the present. In particular, the scheme is 
able to deal with haloes that pop into existence, that grow 
in size by accreting additional background particles, that 
merge with other haloes of comparable size and lose their 
identity, or that fall into a larger halo without being de- 
stroyed completely. The latter case is particularly interest- 
ing, since we expect that a subhalo can survive for some 
time within the larger structure. However, its mass can be 
reduced by tidal stripping, an effect that can eventually com- 
pletely dissolve the structure. 
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4.5 Inclusion of subhaloes in semi-analytic models 

One of the questions we want to address in this study is 
how the inclusion of subhaloes changes the results of semi- 
analytic models. In order to highlight such changes we will 
modify the 'standard' scheme, which is based on the work of 
KCDW, in a minimum fashion when subhaloes are included. 

Note that a large variety of semi-analytic "recipes" are 
conceivable, and that some of the results can depend sen- 
sitively on the specific set of assumptions made, as has 
been recently highlighted by KCDW. Semi-analytic mod- 
els thus cannot eliminate uncertainties arising from poorly 
constrained physics. However, they are ideal tools to ex- 
plore the relative importance of various model ingredients, 
and hence to constrain their relevance for observed trends 
in observational properties. 

We now describe the changes in our semi-analytic 
methodology when subhaloes are included. From here on, 
we will refer to the formalism of KCDW, which only works 
with FOF groups, as the 'standard'-scheme, and to the anal- 
ysis that includes substructure as the 'subhalo'-scheme. We 
note that the word 'standard' is not meant to imply that 
the corresponding procedure has (yet) the status of a well- 
established, widely used method in the field - it is just used 
to refer to the methodology recently developed by KCDW. 

Our most important change concerns the definition of 
the galaxy population at each output time. We define the 
largest subhalo in a FOF-group to host the central galaxy of 
the group, and this galaxy's position is given by the most- 
bound particle in that subhalo. All gas that cools within 
a FOF-group is funneled exclusively to the central galaxy. 
This definition of 'central galaxy' thus corresponds to the 
one adopted in the standard analysis. 

For the population of galaxies orbiting within a halo, 
however, we distinguish between halo-galaxies and satel- 
lite galaxies. Here we have coined yet another term; 'halo- 
galaxies' are attached to the most bound particle of the re- 
maining subhaloes in the FOF-group. These halo-galaxies 
were proper central galaxies in the past, until their halo 
fell into a larger structure. The core of this halo is, how- 
ever, still intact, and thus allows an accurate determination 
of the position of the halo-galaxy within the group. These 
halo-galaxies may still be viewed as 'central galaxies' of their 
respective subhaloes, but they are no longer fed by a cooling 
flow since their subhalo is not the largest within the FOF 
group. 

Finally, when two (or more) subhaloes merge, the halo- 
galaxy of the smaller subhalo becomes a satellite of the rem- 
nant subhalo. These satellites are treated as in the standard 
analysis. Their position is tagged by the most-bound particle 
identified at the last time they were still a halo-galaxy, and 
they are assumed to merge on a dynamical friction timescale 
with the halo-galaxy of the new subhalo they now reside in. 
We need to introduce such satellites in the subhalo-scheme in 
order to account for actual mergers between subhaloes, and 
also because of the finite numerical resolution of our simula- 
tions, which limits our ability to track the orbits of subhaloes 
once their mass has fallen below our resolution limit. It also 
allows us to make direct contact with the standard-scheme 
in the limit of poor resolution. Note that the class of halo- 
galaxies is absent in the standard analysis, where all of these 
galaxies are treated as satellites. 



Table 2. Numerical parameters adopted for our semi-analytic 
models, a is the star formation efficiency, e the efficiency of feed- 
back by supernovae, and /t, the baryon fraction. 



Model 


a 


e 


fb 


ejection 


0.05 


0.10 


0.15 


retention 


0.05 


0.15 


0.15 



In the subhalo scheme, we define the virial mass M V1I of 
a subhalo simply as the total mass of its particles. For the 
background subhalo, we then define virial radius and virial 
velocity by assuming that the halo has an overdensity 200 
with respect to the critical density. For other subhaloes we 
keep the virial velocity and the halo's dynamical time at the 
values they had just before infall. 



5 RESULTS 

5.1 Tully-Fisher relation 

We use the velocity-based /-band Tully-Fisher relation 

Mi - 5 log h = -21.00 - 7.68 (log W - 2.5) (14) 

of Giovanelli et al. (1997), and the requirement of a gas mass 
of ~ 8 x 10 9 /i _1 M Q in 'Milky- Way' haloes, to normalize our 
models. We consider two variants for the implementation 
of feedback, the 'ejection' model, where gas is blown out of 
small haloes, and the 'retention' model, where reheated gas 
is always kept within the halo. 

In Figure 6, we show the best-fit Tully-Fisher relations 
obtained for these two models, applied to the S2-cluster us- 
ing the 'subhalo' and the 'standard'-schemes. In the plots, 
we only considered central galaxies of haloes that are pe- 
ripheral to the cluster, but that are not contaminated by 
heavier boundary particles. We also applied a morphologi- 
cal cut, 1.2 < Mbuigc — Mtotai < 2.5, approximately selecting 
Sb/Sc galaxies. In Table 2, we list the model parameters thus 
obtained. In the following, we will use the same set of pa- 
rameter values also for our other cluster simulations, and for 
the 'standard' semi-analytic scheme. 

The Tully-Fisher relations we obtain exhibit remarkably 
small scatter. This is partly a result of the tight coupling we 
assumed between the sizes and circular velocities of the disks 
of spiral galaxies and the masses of their dark haloes. Note 
that additional scatter can be expected from the distribu- 
tion of spin parameters of dark haloes, which gives rise to 
variations of the disk sizes associated with a halo of a given 
mass (Mo et al. 1998). 

The ejection model fits the slope of the observed TF- 
relation relatively easily. However, the retention model is less 
effective in suppressing star formation in low mass haloes. 
For the same value of e, the retention scheme therefore pro- 
duces a shallower Tully-Fisher relation than the ejection 
model. As a result, a larger value of e is needed to bring 
the retention model in agreement with the observed steep- 
ness of the TF-relation. This strong feedback reduces the 
overall brightness somewhat, an effect that could be eas- 
ily compensated for by slightly larger values for a or fb. 
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Figure 6. The /-band Tully-Fisher relation for Sb/Sc galaxies in the S2 simulation. The two panels on top have been obtained with the 
subhalo-schcmc using ejection and retention feedback, respectively, while the two panels on the bottom show the equivalent plots using 
the standard-scheme of KCDW. The galaxies have been selected as central galaxies of uncontaminated haloes in the periphery of the 
cluster. The solid line represents the recent measurement by Giovanelli et al. (1997). 



Note that a feedback efficiency of e = 0.1 means that feed- 
back will be quite efficient in haloes of virial velocity below 
Vsn = 183kms~ 1 (e/0.1) 1//2 . In such haloes, one will have 
AM ro hcat > AM*, i.e. the mass of reheated gas exceeds that 
of newly formed stars. 

It is also interesting to compare the Tully-Fisher rela- 
tions obtained for the different cluster simulations (see Fig- 
ure 7) . For the sake of brevity, we restrict the comparison to 
the subhalo scheme with ejection feedback. In general, there 
is good agreement between the two variants of semi-analytic 
modeling, and between the four different simulations, de- 
spite their large differences in numerical resolution. For the 
same choice of free parameters, the slopes and zero-points 
of the Tully-Fisher relations agree quite well. There may 
be a weak trend towards fainter zero-points in the sequence 
SI to S4. Note that the scatter in the Tully-Fisher rela- 
tion increases at low velocity widths because feedback has 
a stronger effect on galaxies with low V c . A similar trend 
in the observed scatter has been found by Matthews et al. 
(1998). 



5.2 Cluster luminosity function 

In Figure 8, we show the B-band luminosity function of 
the S2-cluster obtained with the new 'subhalo' methodol- 
ogy, and we compare it to the result of the 'standard' semi- 
analytic recipe of KCDW. We plot the number of cluster- 



galaxies in bins of size 0.5 mag, and we fit Schechter func- 
tions to the counts. The standard prescription results in 
a relatively steep slope of —1.31 at the faint-end, and the 
"knee" of the Schechter function is not well defined. This is 
just another reincarnation of a well-known problem in many 
previous semi-analytic studies, which predicted too many 
galaxies both at the faint and the bright end of the the 
luminosity function, i.e. the shape was not curved enough. 
These deficiencies can be partly cured by invoking additional 
physics like dust obscuration, or by employing models with 
very strong feedback processes. However, most semi-analytic 
studies have not been successful in simultaneously achieving 
good agreement with the Tully-Fisher relation and the lu- 
minosity function. 

Compared to the field, the luminosity functions of clus- 
ters tend to be considerably steeper, and a slope of —1.31 
can be accommodated with existing data. However, the 
standard-scheme also produces a brightest cluster galaxy 
with luminosity in excess of most normal cD galaxies. For 
example, in the S2 cluster, the standard-scheme produces a 
central galaxy with B-magnitude —23.9. 

In contrast, the luminosity function obtained with the 
subhalo formalism has a more curved shape. There are 
more L*-galaxies, resulting in a flatter faint-end slope of 
a = —1.21, and in a more pronounced "knee" , which is much 
better fit with a Schechter function. In addition, the lumi- 
nosity of the central galaxy is reduced. The overall shape 
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Figure 7. The four panels compare the /-band Tully-Fisher relations obtained for the SI, S2, S3, and S4 simulations using the subhalo- 
scheme and ejection feedback. The models with retention feedback, and the models using the 'standard' scheme, are able to match the 
observed Tully-Fisher relation about equally well. 



of the resulting luminosity function is in reasonably good 
agreement with observed cluster luminosity functions, as 
shown by a comparison with the composite luminosity func- 
tion of Trentham (1998), which is indicated with symbols in 
Figure 8. 

What causes this difference between the subhalo-scheme 
and the standard scheme? Note that the excessive brightness 
of the central galaxy in the latter model is unlikely to be 
caused by an overcooling problem. We have already cut-off 
the cooling flow for central galaxies in haloes with a virial 
velocity larger than 350 km s -1 . Furthermore, the cooling 
model was essentially the same for both schemes, which is 
also reflected in their similar overall mass-to-light ratios. 

The most important difference between the 'standard' 
semi-analytic scheme and the 'subhalo'-methodology is the 
explicit tracking of subhaloes, which allows more faithful 
modeling of the actual merging rate in any given halo. Recall 
that one critical assumption in the study of KCDW is that 
the time of survival of a satellite can be estimated using a 
simple dynamical-friction formula. However, this description 
is crude and the excessive brightness of the central galaxies 
may result from an overestimate of the overall merging rate, 
or from merging the 'wrong' galaxies. Recall that a second 
important assumption of KCDW has been that the position 
of a satellite galaxy can be traced by a single particle iden- 
tified as the most-bound particle of the satellite's halo just 
before it was accreted by a larger system. 

Before we investigate these assumptions in more detail, 



we plot in Figure 9 the cluster luminosity functions for the 
SI, S2, S3 and S4 clusters, obtained with the subhalo scheme 
and ejection feedback. All four simulations produce luminos- 
ity functions which can be well described by Schechter func- 
tions with a well defined cut-off. Their shape is in fact quite 
close to the obervational result by Trentham (1998), shown 
with symbols in each of the four panels. The vertical nor- 
malization of his composite observational result is arbitrary, 
and we have set it to match the total luminosity of our clus- 
ter^). It is however the same in all four panels. Agreement 
of the luminosity functions between the four simulations is 
quite good, with the simulations of higher mass resolution 
able to probe the luminosity function to ever fainter mag- 
nitudes. It is interesting to note that as the resolution of 
the simulations increases, a larger and larger fraction of the 
galaxies have their own dark matter halo (this is indicated 
by a dashed line in the plot). In S4, almost all bright galax- 
ies still have an associated dark matter substructure, and 
it seems likely that in still larger numerical simulations the 
correspondence would be perfect. 

However, when comparing the luminosity functions 
with each other, there is a trend of decreasing brightness 
of the first ranked cluster galaxies with increasing resolu- 
tion. We think that this is a reflection of inaccuracies in 
estimates of merging timescales, which become more serious 
in simulations with lower resolution. This effect may also 
be responsible for the weak trend in the zero-points of the 
Tully-Fisher relations (Figure 7). 
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Figure 8. The histograms show the B-band cluster luminosity function obtained for the S2 simulation using the subhalo scheme (left) 
and the standard scheme (right). We here show results for the ejection feedback, and note that the results for the retention model are 
very similar. The symbols and error bars give the composite luminosity function found by Trentham (1998) as a weighted mean of 9 
clusters, all having luminosity functions that are individually consistent with the composite one. Note that we adjusted the vertical 
normalization of the composite function (which is arbitrary) to obtain a good fit to the total luminosity of the cluster. The solid lines are 
Schechter function fits to the histograms of bin-size 0.5 mag. For the subhalo model, the faint end slope is a = —1.21, and the turn-off 
is at Mi, ~ —21.8. The standard scheme results in a steeper slope of a = —1.31, and the characteristic magnitude is not well defined. 
Note that the standard model produces a central galaxy of magnitude -23.9, which appears too bright even comparing to the brightest 
cD galaxies. 



We now test whether this problem is responsible for 
the difference between the subhalo scheme and the standard 
formalism. To this end we try to answer the question: How 
many of the galaxies present in subhaloes of S2 at z = arc 
prematurely merged with the center in the standard scheme? 

To do this, we follow each subhalo of the cluster back in 
time along its merging history until it is a main halo itself 
for the first time. The corresponding FOF-halo will host a 
central galaxy in the standard scheme, and the descendant 
of this galaxy at z = should directly correspond to the 
galaxy of the originally selected subhalo. 

Among the 494 subhaloes identified in the S2-cluster 
at z = 0, we find in this way that 39 of them do not have 
a directly corresponding satellite galaxy in the 'standard'- 
formalism. The satellites that should correspond to these 
39 subhaloes have prematurely disappeared by merging 
processes with the central galaxy, because their merging 
timescales have been underestimated. Note that the total 
number of mergers with the central galaxy in the standard 
scheme is 116, while this number is 132 in the subhalo for- 
malism. This suggests that the overall rate of mergers is 
not too high in the standard scheme. However, in the stan- 
dard scheme the central galaxy accretes 24 galaxies with 
stellar masses of more than 10 10 /i _1 Mq, among them 8 
galaxies with stellar masses larger than 10 11 ft _1 M0. On the 
other hand, in the subhalo formalism there is only 1 merger 
involving more than 10 11 /i _1 Mq, and 11 with more than 
10 10 /i _1 Mq . This means that a larger number of very bright 
galaxies is merged with the central galaxy in the standard 
scheme. We also note that the galaxies in the 39 subhaloes 
that appear to have been merged prematurely in the stan- 



dard scheme tend to be quite bright. If we take the results 
for the subhalo-model and add the luminosities of the cor- 
responding halo-galaxies to that of the central galaxy, its 
brightness increases from -21.9 to -24.1 mag, quite close to 
the -23.9 obtained in the standard-scheme. It thus appears 
that the excessive brightness of the central cluster galaxy in 
the standard scheme is mainly caused by an underestimate 
of the merging timescales for some fraction of the bright 
galaxies. 

We have also tested how well the position of subhaloes 
corresponds to those of single particles used to track satel- 
lites in the standard scheme. We find that a fraction of 86.8% 
of the subhaloes in S2 still contain the most-bound particle 
that is used in the standard-scheme to track the position of 
the satellite. This number is 86.6% in SI, 82.2% in the S3 
cluster, and 76.3% in the S4 cluster. Note that these num- 
bers have been computed for all the subhalos in each of the 
simulations. To the mass resolution limit of the SI cluster, 
the corresponding success rates are 93.6%, 89.4% and 96.8% 
in S2, S3, and S4, respectively. Using just a single particle 
identified at a time before a structure was accreted onto the 
cluster can thus provide a good estimate of the subhalo's 
position within the cluster halo at later times. 

It thus seems clear that the difference between the lu- 
minosity functions of the standard and subhalo schemes 
is largely caused by differences in the assumed merger 
timescales of galaxies that fall into the cluster. Recall that 
for the choices we adopted in equation (9), the dynamical 
friction formula is essentially the current dynamical time of 
virialized haloes times the mass ratio of the halo and its in- 
falling satellite. Large estimates of merger timescales arise 
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Figure 9. The B-band cluster luminosity functions obtained for the SI, S2, S3, and S4 simulations using the subhalo scheme and ejection 
feedback. The solid lines are three-parameter Schechter function fits to the histograms of bin-size 0.5 mag. The formal values for faint-end 
slope and cut-off of these fits arc: SI (a = -1.22, M* = -22.8), S2 (a = -1.21, = -21.8), S3 (a = -1.21, M* = -21.7), and S4 
(a = —1.24, M* = —22.0). Symbols with error bars give the composite observational result by Trentham (1998). Dashed histograms give 
the luminosity functions of subhalo galaxies alone, i.e. galaxies that still have an associated dark matter subhalo. 



when the mass of the satellite is much smaller than that of 
the halo. In fact, small satellites will typically have merger 
timescales that are much larger than a Hubble time, while 
large infalling haloes are predicted to merge with the cen- 
tral galaxy quickly. As a result, the standard scheme appears 
to merge some of the associated bright galaxies 'too early', 
making the central galaxy excessively bright. In the sub- 
halo scheme, these galaxies are still around and populate 
the 'knee' of the luminosity function. Note however that the 
dynamical friction estimate in the simple form applied here 
does not include the effects of tidal truncation and disrup- 
tion, which presumably act to reduce the lifetime of small 
satellites, but may increase the lifetime of infalling massive 



satellites. In fact, the results of Tormen et al. (1998) indicate 
that the dependence of the merger timescale on the mass ra- 
tio, measured just before the satellite falls in, is weaker than 
linear. 

Presently it is unclear whether an improved parameter- 
ization of the merger timescales in the standard scheme can 
produce results as good as those obtained by explicitly fol- 
lowing the subhaloes. However, a more detailed investigation 
of the probability distribution of actual merging timescales 
may well make this possible, and offers the prospect of im- 
plementing more accurate satellite merging in the 'standard' 
semi-analytic scheme, even when simulations with relatively 
poor resolution are used. 
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Figure 10. Cumulative number of galaxies in the cluster as a 
function of radial distance to the cluster center. In the top panel, 
we count all galaxies with M^, < — 16, binned into three classes 
of different Hubble type. In the bottom panel, only the bright 
galaxies with M-g < — 18 are shown. In both panels, the shape of 
the cumulative dark matter mass distribution is indicated by a 
thin solid line, for comparison. The results shown are for the S4 
cluster using ejection feedback. 



5.3 Morphology density relation 

Dressier (1980) has shown that the relative frequency of el- 
liptical galaxies is higher in denser environments. In partic- 
ular, in the cores of clusters, spirals are quite rare, while 
they are the dominant type in the field, and in the Universe 
as a whole. Whitmore et al. (1993) have argued that this 
morphology-density relation reflects a more fundamental 
morphology-clustercentric relation; the correlation between 
morphology and clustercentric separation seems tighter than 
that between morphology and projected density. This asser- 
tion is still controversial. 

In Figure 10, we show the cumulative number of galax- 
ies of different Hubble type as a function of distance from 
the cluster center. The top panel gives all the galaxies with 
Mb < —16, while in the bottom panel we just show the 
bright galaxies with Mb < —18. From the cumulative distri- 
bution it can be inferred that the three-dimensional density 



profile of the ellipticals has a steeper slope than that of the 
spirals, i.e. they are more concentrated towards the center. 
Note that the thin solid line gives the shape of the cumula- 
tive dark matter profile for comparison. It is also interesting 
to note that the bright galaxies in the center are primarily 
ellipticals, while their mean Hubble type shifts to later types 
towards the outskirts of the cluster. 

These trends in the morphological mix of galaxies can 
be more clearly seen in Figure 11, where we show the 
morphology-clustercentric relations obtained for the S2 and 
S4-clusters, using the subhalo-scheme with ejection feed- 
back. Note that the fraction of elliptical galaxies strongly 
rises towards the center of the cluster, while that of the spi- 
rals declines accordingly. The quantitative strength of these 
trends is in reasonable agreement with the results of Whit- 
more et al. (1993), although there are small differences in 
detail. For the classification of galaxies, we applied cuts in T- 
space that count galaxies with < T < 5 as SO's, and lower 
(higher) types as ellipticals (spirals), respectively. The re- 
sulting relative populations of galaxies of different types are 
consistent with the observed morphological mix, and largely 
independent of the magnitude cut-off for the sample. Note 
however that in the standard observational classification a 
slightly different cut in T-space is used for SO's. This may 
just reflect the difficulty to unambiguously define the tran- 
sitional type SO in our coarse morphological classification 
scheme, or it may mean that additional physical processes 
are at work to create real SO galaxies. 

Recall that our morphological modeling is solely based 
on the merging history of galaxies. Ellipticals are formed in 
major mergers, which occur more frequently in higher den- 
sity environments. This simple model suffices to establish a 
pronounced morphology density relation. This shows that 
the morphology density relation is built-in at a very funda- 
mental level in hierarchical theories of galaxy formation. 

In passing we note that in the standard recipe, where 
satellites are just traced by single particles once they have 
fallen into a larger halo, the morphology-clustercentric re- 
lation is also present, although it is not quite so well de- 
fined. The stronger clustering of early type galaxies was also 
found by KCDW, Diaferio et al. (1999) and more recently 
by Diaferio et al. (2000). Using techniqes somewhat sim- 
ilar to ours, Okamoto & Nagashima (2000) also found a 
morphology-density relation. All these papers argued that 
additional processes are required to understand the SO pop- 
ulation. 

Figure 11 also shows that the fraction of SO's in our clus- 
ter shows a weaker dependence on cluster-centric distance 
than the observations. This also suggests that SO's experi- 
ence additional processes in clusters that are not included in 
our analysis. For example, disk galaxies orbiting in a clus- 
ter experience high-speed encounters with other galaxies. 
Together with global cluster tides this can drive a morpho- 
logical evolution towards spheroids, a process termed galaxy 
harassment (Moore et al. 1996, 1998). It has been suggested 
that the spiral galaxies seen abundantly in clusters at mod- 
erate redshift are harassed and slowly transformed to SO's 
at the present time. Such an effect might explain our deficit 
of SO galaxies inside the cluster. 
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Figure 11. Morphological mix of galaxies as a function of clustercentric radius. From top to bottom, the symbols in the three panels show 
the relative fraction of elliptical, SO, and spiral galaxies in spherical shells around the cluster center, with error bars given by counting 
statistics. Filled symbols are for galaxies brighter than Mb = —17, and hollow symbols for a sample selected at Mb = —18.5. Filled 
circles show the observational results of Whitmore et al. (1993). Note that we here follow these authors in letting the radius decrease 
to the right, i.e. the cluster center is found on the right side of the diagrams. The results shown here are for the S2 and S4 simulations, 
using the subhalo modeling with ejection feedback. Results for the other two simulations are similar. 



5.4 Cluster mass-to-light ratio 

The observed mass-to-light ratios of clusters of galaxies are 
known to be much larger than those of individual galaxies, 
and this fact has been recognized early on as strong evidence 
for the existence of large amounts of dark matter in clusters. 
Typical measured values for the cluster mass-to-light range 
from T s = 200 ft T© to T s = 400 h Tq in the B-band. For 
the Coma cluster, Kent & Gunn (1982) measured Ts = 
360 ft T©, while X-ray data seem to point to a lower value. 
For example, Cowie et al. (1987) inferred a mass-to- visual- 
light ratio for Coma as low as Ty = 180 ft Y©. 



The galaxy population constructed for the S4 cluster 
using the subhalo model has a total magnitude of Mb = 
-26.04. For its total mass of 8.36 x 10 14 ft _1 M , the cluster 
mass-to- light ratio is thus Ts = 420 ft T©. In the F-band, 
we obtain My = —26.97 and Ty = 324 ft T©. These mass- 
to-light ratios fall on the high side compared to the mean of 
observational results, although there are also measurements 
that are as large, or larger, as our values. For example, Kent 
& Gunn (1983) obtained Ty = 600 ft T© for the Perseus 
cluster. 

Small changes of the model parameters can, however, 
enhance the overall brightness of the cluster, and thus re- 
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Figure 12. B-band mass-to-light ratio of halo galaxies as a func- 
tion of distance to the cluster center (small symbols), expressed 
in units of the mass-to-light ratio of the cluster as a whole. For 
the masses of the galaxies we here simply took the dark matter 
masses of their corresponding subhalocs. The large symbols give 
the median mass-to-light ratios in radial bins around the clus- 
ter center. There is a clear trend of decreasing M/L towards the 
cluster center, largely reflecting the mass loss of subhaloes due to 
tidal truncation. Note however that early type galaxies appear to 
have a systematically lower mass-to-light ratio than late types. 



duce the mass-to-light ratios to a desired value. One simple 
possibility is to change the conversion factor between the 
virial velocity and the circular velocity. If the circular veloc- 
ity of a disk in a dark matter halo of given virial velocity is 
on average 15% larger than we assumed, the normalization 
to the Giovanelli et al. (1997) Tully-Fisher relation results 
in an overall increase in brightness of 0.5 mag, and a corre- 
sponding reduction of the mass-to-light ratio to 63% of its 
old value. 

Dust obscuration can account for effects of similar size. 
It is quite likely that the observed average brightness of 
galaxies in the Giovanelli et al. (1997) sample is somewhat 
reduced by dust, even though the /-band is less affected by 
dust than bluer wavelengths. Since we do not correct for dust 
extinction in this work, our fit to the observed Tully-Fisher 
relation is expected to underestimate the stellar masses of 
these spirals. Correcting this would produce larger stellar 
masses for all galaxies at the present epoch and consequently 
lower mass-to- light ratios for the cluster. 

We may also ask whether we detect systematic varia- 
tions of the mass-to-light ratio of subhalo galaxies as a func- 
tion of their position in the cluster. In Figure 12 we show the 
B-band mass-to-light ratio of halo galaxies as a function of 
radial distance to the cluster center. Here we took the dark 
matter masses of the corresponding subhaloes to compute 
M/L ratios. The clear trend of decreasing M/L towards the 
cluster center is caused primarily by mass loss through tidal 
truncation. However, for a given distance from the cluster 
center, there is also a systematic variation of mass-to-light 
ratio with morphological type. At each luminosity, early type 
galaxies tend to have lower mass halos than late types, re- 
flecting their earlier incorporation and more effective strip- 
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Figure 13. Faber Jackson relations for cluster ellipticals (cir- 
cles) and field/group ellipticals (filled stars). We here plot the 
f?-band luminosity of subhalo-galaxies versus an estimate of the 
one-dimensional stellar velocity dispersion. The latter was taken 
to be 0.9/\/3 times the measured 3D- velocity dispersion of the 
galaxies' dark matter subhalo identified in the simulation. The 
factor l/vo converts from 3D to line-of-sight velocity dispersion, 
and the factor 0.9 is needed to bring the zeropoint in agreement 
with the observed correlation. The obervational Faber-Jackson 
relation is shown as a solid line, and is adopted from the fit of 
Bender et al. (1996) to data for Coma and Virgo ellipticals. 

ping within the cluster. All the galaxies have mass-to-light 
ratios substantially smaller than the cluster as a whole. 

5.5 Faber-Jackson relation 

For elliptical galaxies, there is a well-known scaling rela- 
tion between luminosity and central velocity dispersion, the 
so-called Faber-Jackson relation. In our formalism, we ex- 
plicitly follow dark matter subhaloes that orbit inside larger 
group- or cluster-sized haloes. While the mass of these sub- 
haloes decreases strongly over time by tidal truncation, their 
velocity dispersion remains relatively stable until they finally 
disrupt. Thus the velocity dispersion of a heavily truncated 
subhalo still reflects that of the original halo just before it 
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fell into the cluster. It is plausible that the stellar veloc- 
ity dispersions of elliptical galaxies correlate with the dark 
matter dispersions of their haloes (Rix et al. 1997). We may 
thus hope to find a Faber- Jackson-like relation if we plot the 
luminosity of our elliptical galaxies aginst the dark-matter 
velocity dispersion of their surrounding (sub)haloes. 

In Figure 13 we show the resulting _B-band Faber- 
Jackson relation for the S2 and S4 simulations. We here 
converted the measured 3D velocity dispersions of the dark 
matter of the subhalos to ID dispersions, and we multiplied 
the resulting values by a factor 0.9 to obtain an estimated 
stellar velocity dispersion. It is remarkable that the zero- 
point, slope and scatter are consistent with the observed 
relation for cluster ellipticals, which is here adopted from 
Bender et al. (1996) for the Virgo and Coma clusters and 
drawn as a solid line. It is reassuring that the kinematics of 
the dark matter substructure has indeed the right properties 
to explain the Faber- Jackson relation. Note that some of the 
bright elliptical galaxies in this plot have heavily stripped 
subhaloes of relatively small mass, but the dark matter ve- 
locity dispersion of the residual subhalo is still large enough 
to put them at about the right place according to the Faber- 
Jackson relation. 

5.6 Luminosity segregation 

We have already seen that inside the cluster the number dis- 
tribution of luminous elliptical and spiral galaxies does not 
follow the mass distribution of the dark matter. Does this 
also hold for the total light emitted by the cluster galaxies? 
In Figure 14, we show the cumulative luminosity profile of 
the cluster, and compare it to the dark matter mass profile. 
The luminosity profile is shallower than that of the mass, 
showing that the light is actually more concentrated than 
the mass, or phrased differently, the cumulative mass-to- 
light ratio decreases towards the cluster center. Note that 
the 'steps' in the light profile are caused by the finite num- 
ber of galaxies and our assumption that they emit as point 
sources. 

A related question is that of luminosity segregation, 
i.e. are more luminous galaxies more strongly concentrated 
than less luminous ones? Note that clustering strength is a 
strong function of circular velocity in CDM models (White 
et al. 1987), so some form of luminosity segregation should 
perhaps be expected. Since early-type galaxies tend to be 
brighter than late-type galaxies, luminosity segregation may 
also be seen as a consequence of the morphology density re- 
lation. However, it is then not so clear which of these two 
correlations is more fundamental. In fact, they might have 
a common origin in the local dynamics of the cluster core 
with its frequent mergers. 

Observational evidence for luminosity segregation is 
still relatively sparse and controversial, primarily because 
of the difficulty in obtaining sufficiently large samples of 
faint cluster galaxies. Kashikawa et al. (1998) have recently 
studied the dependence of luminosity segregation on mor- 
phological type in the Coma cluster. In order to disentan- 
gle luminosity segregation from the morphology-density re- 
lations it is critical to understand whether it occurs within 
a given morphological type. Kashikawa et al. (1998) found 
that galaxies with high central concentration (early types) 
show signs of luminosity segregation, while galaxies of low 
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Figure 14. Luminosity segregation in the cluster S4. We here 
show the cumulative luminosity profile of B-band light of cluster 
galaxies (thick solid lines), and compare it to the cumulative dark 
matter mass profile (thin sold lines). The latter has been scaled 
by the overall mass-to-light ratio of the cluster to put the curves 
in the three panels on a common scale. While the top panel shows 
all the galaxies, we have devided them by morphological type into 
two equally luminous groups in the middle and bottom panels. In 
each panel, we further split the galaxies into a luminous (dashed 
lines) and a faint (dotted lines) group such that each of the groups 
produces just half of the total luminosity in the panel. Finally, we 
shifted the total luminosity profile in each panel by a factor of two 
to align the curves at the virial radius. 



central concentration (late types) do not exhibit any strong 
segregation. 

In the two bottom panels of Figure 14, we investigate 
this issue for our cluster galaxies. We have divided the galax- 
ies into early and late types, and we split each of the groups 
into a high and a low luminosity sample, such that the total 
luminosity in each group was divided into two equal parts. 
Interestingly, only the light of the bright elliptical galaxies 
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Figure 15. One-point velocity dispersion of galaxies (filled circles) brighter than Mb = —17 as a function of distance to the cluster 
center. The velocity dispersion of the dark matter is given by the solid lines in the four panels. The error bars mark the 68 per cent 
confidence intervals based on counting statistics. We have further split up the galaxies into a blue and a red sample of equal size, based 
on their B — V color. The blue galaxies are shown as triangles, the red ones as boxes. Note that there are no very blue galaxies near the 
cluster center, hence triangular symbols are absent at small radii. 



is more concentrated than the mass of the cluster, while the 
low luminosity early types, and the spiral galaxies are dis- 
tributed like the mass. We thus find evidence for luminosity 
segregation of early type galaxies. This is in agreement with 
the findings of Kashikawa et al. (1998). Notice however that 
Diaferio et al. (2000) found no clear evidence of luminosity 
segregation in their analysis of the GIF simulations. 



5.7 Velocity dispersion of galaxies 

One common technique to estimate masses of real clusters 
of galaxies is based on the virial theorem, requiring mea- 
surements of the line-of-sight velocity dispersion and of the 
projected separations of galaxies (Heisler et al. 1985). How- 
ever, the accuracy of this technique can be compromised if 
cluster galaxies are biased tracers of the dark matter veloc- 
ity distribution. Such systematic differences in the velocity 
fields will usually give rise to errors in the mass estimates 
of clusters. In fact, there have been controversial claims in 



the literature whether such a bias exists, and whether it is 
positive or negative. Since there is currently no consensus 
on this issue, we briefly address it here. In the following, we 
will focus on the one-point velocity bias, defined as the ratio 
b v = Cgai/cDM °f g aiax y and dark matter velocity disper- 
sions. We will refer to values of b v greater (less) than one as 
positive (negative) velocity bias. 

In the early simulation work of Carlberg (1994) and 
Frenk et al. (1996), a negative velocity bias of up to 20- 
30% was reported. However, using the 'standard' scheme 
of semi-analytic modeling applied to the GIF-simulations, 
Diaferio et al. (1999, 2000) found a positive velocity bias in 
the outskirts of clusters, a result which can be explained as 
being due to the recent infall of star-forming blue galaxies. 

Recent numerical studies with sufficient resolution to 
resolve substructure have led to somewhat conflicting re- 
sults. Colin et al. (2000) found a substantial positive velocity 
bias, except perhaps in the innermost region of their clus- 
ter. A similar result, albeit of weaker strength, was obtained 
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by Okamoto & Habe (1999). On the other hand, Ghigna 
et al. (1998) did not find any significant velocity bias in 
their analysis of cluster subhaloes. Similarly, Klypin et al. 

(1999) found no velocity bias on cluster mass scales, but a 
significant anti-bias for small groups of galaxies. In a more 
detailed analysis of this issue, Ghigna et al. (1999) recently 
reported a small positive velocity bias for subhaloes in the 
innermost region of their clusters, although the overall sig- 
nal was considerably smaller than that found by Colin et al. 

(2000) . 

In light of this theoretical uncertainty, it is interest- 
ing to examine the velocity bias of our cluster simulations. 
In Figure 15, we plot the velocity dispersion of our model 
galaxies in the clusters S1-S4 as a function of radius. We 
used logarithmic bins and computed 68%-confidence inter- 
vals under the assumption that the galaxies in each bin are 
drawn from an underlying Gaussian velocity distribution. In 
this case, the estimated dispersions are expected to follow 
Xn-i distributions, where n is the number of galaxies in each 
bin. 

Upon comparison with the dark matter velocity dis- 
persion (solid lines) it is seen that we detect evidence for 
a small negative velocity bias in the central region of the 
cluster. The significance of this result is small for a single 
bin, but adjacent bins are independent here, and 'noise' be- 
tween the four realizations is also uncorrelated, giving rise 
to a quite significant finding when the individual results are 
combined. We thus conclude that our models show evidence 
for a negative velocity bias in the innermost ~ 300/i _1 kpc 
of the cluster, while further outside, velocity bias appears 
to be small, if present at all. Such a lower velocity disper- 
sion in the center appears to be consistent with a number of 
observational studies (e.g. Biviano et al. 1992; Adami et al. 
1998) that have found systematically smaller velocities for 
the brightest galaxies in the cores of clusters, although the 
observational evidence for this has not so far been fully con- 
vincing. 

In Figure 15, we have also split up the galaxies into a 
blue and a red sample according to their B — V color. Inter- 
estingly, it is seen that the blue galaxy sample shows larger 
velocity dispersions in the radial range 300 - lOOO/i^kpc. 
This is in agreement with the findings of Diaferio et al. (1999, 
2000). Since these galaxies experienced relatively recent star 
formation, they have mostly just fallen into the cluster. As 
a result their orbits have systematically larger apocentres 
than those of red galaxies in the same region of the clus- 
ter. Note that if one selects the most massive subhaloes, one 
also obtains a population that has just fallen into the cluster 
and is biased towards its outskirts. This may explain why 
studies focusing on subhaloes have reported signs of positive 
velocity bias (Colin et al. 2000). 

5.8 B — V color distribution 

In Figure 16, we show the B — V color distribution of galax- 
ies brighter than —19.7 mag in the S2 simulation. A cor- 
responding plot has been shown by KCDW, and we here 
obtain similar result. The distribution of B — V colors of 
field galaxies is bimodal, with two peaks at B — V ~ 0.7 
and B — V ~ 1.0. By plotting the color distribution for in- 
dividual morphological types, it becomes apparent that this 
dichotomy mainly reflects the differences in star formation 
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Figure 16. The B — V color distribution of galaxies brighter 
than Mb = —19.7 mag in the 'field' region of the S4 cluster 
simulation (shaded histogram). The dichotomy of the distribu- 
tion mainly reflects the difference in the stellar populations of 
elliptical and spiral galaxies. Spirals are significantly bluer due to 
the their recent star formation, while ellipticals have older stel- 
lar populations, resulting in redder colors. The thick histogram 
gives the observed color distribution as determined by KCDW 
from the RC3 catalogue (de Vaucoulers et al. 1991). Finally, the 
dashed line outlines the B — V color distribution of galaxies in the 
cluster. Clearly, these galaxies tend to be much redder than the 
field sample. Even most of the cluster-galaxies that are classified 
as spirals are dominated by relatively old stellar populations. 



history between elliptical and spiral galaxies. The recent star 
formation in spirals makes them blue, while the older stellar 
populations of ellipticals give them redder colors. The color 
distribution of field galaxies is in good agreement with the 
RC3 sample (de Vaucoulers et al. 1991) of local galaxies. In 
Figure 16 we also plot the color distribution of the galax- 
ies in the cluster (dashed histograms). These galaxies tend 
to be much redder, reflecting their substantially older stel- 
lar populations. Note that since we have not included dust 
in our models, the apparent agreement in Figure 16 proba- 
bly means that the populations of late-type galaxies in our 
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model are systematically older than those of the correspond- 
ing galaxies in the RC3. 



6 DISCUSSION 

In this study, we have used cosmological N-body simula- 
tions combined with semi-analytic techniques to construct 
the galaxy population of a rich cluster of galaxies. The very 
high resolution of our simulations allowed us to extend the 
methodology for galaxy formation of KCDW to the regime 
of substructure within virialized systems. Our goal in this 
work has therefore been twofold. We wanted to develop the 
necessary technical machinery for a galaxy formation scheme 
that works with subhaloes, and we wanted to compare its 
results with those obtained with the 'standard' techniques 
by KCDW. 

The detection of subhaloes within haloes is a technically 
difficult problem, and has only been addressed very recently. 
Several working algorithms have been described in the lit- 
erature, but none seems ideal. In this paper, we presented 
our new subhalo finding algorithm, SUBFIND. It relies on 
the particle positions and velocities at a single output time, 
and it reliably identifies locally overdense, self-bound parti- 
cle groups within larger systems. SUBFIND can also detect 
hierarchies of 'haloes within haloes', a feature that is not 
possible with alternative techniques. We find examples of 
such a hierarchy in our highest resolution simulation. 

The mass spectrum of subhaloes in our cluster appears 
to be close to a power-law. Our set of simulations allowed a 
direct study of resolution effects, and it is interesting to note 
that the lower-resolution simulations predict the right abun- 
dance of subhaloes for any given mass above their respective 
resolution limit. 

We have shown how the detected subhaloes can be 
traced from simulation output to output. Just like KCDW, 
we analysed 51 simulation snapshots, logarithmically spaced 
in expansion factor from z — 20 to z — 0. This large num- 
ber of output times together with the very high resolution 
of our simulations allowed a measurement of the merging 
history of the dark matter in unprecedented detail. For ex- 
ample, we detected almost 4700 subhaloes in the final halo 
of the S4 cluster, and we found hundreds of mergers between 
subhaloes orbiting inside the progenitor halo of the cluster, 
even though the rate of genuine subhalo-subhalo mergers is 
relatively low. 

Using a small set of modifications, we have adapted 
the semi-analytic scheme of KCDW to include subhaloes 
in the analysis. In both schemes, the agreement with the 
observed Tully-Fisher relation is very good. However, the 
inclusion of subhaloes results in a substantial improvement 
of the cluster luminosity function of the models. In the stan- 
dard scheme, the first ranked cluster galaxies usually turn 
out to be too bright, while this problem does not occur in 
the subhalo-scheme. We have shown that this is mainly due 
to inaccuracies in the estimated merging timescales in the 
standard scheme, where too many bright galaxies are prema- 
turely merged with the central galaxy. The direct tracing of 
subhaloes until their eventual disappearance allows a more 
faithful estimate of the actual merger rate within haloes. As 
a result, the luminosity function becomes more curved, and 
develops a well defined knee and a flatter faint-end slope. 



In fact, in the highest resolution simulations we have car- 
ried out, essentially all the bright galaxies in the cluster can 
still be followed in terms of their individual dark matter 
subhaloes. 

The subhalo-scheme also provides accurate positions 
and velocities for galaxies orbiting in the cluster halo. We 
have shown that our simple morphological modeling gives 
rise to a morphology-clustercentric relation that is quali- 
tatively in good agreement with observations. Towards the 
center of the cluster, the morphological mix of galaxies be- 
comes gradually dominated by ellipticals, while the contri- 
bution of spirals strongly declines. Note that the morphol- 
ogy of our model galaxies is primarily determined by their 
merging history. A morphology density relation arises quite 
naturally in hierarchical theories of galaxy formation. We 
have also found that red galaxies near the cluster center can 
be expected to have a negative velocity bias, i.e. they move 
more slowly than nearby dark matter particles. For infalling 
blue galaxies we detect a small positive velocity bias. 

In summary, we find that our subhalo-model for the 
galaxy population of the cluster produces results that are 
in good agreement with a variety of data. The cluster- 
luminosity function has a reasonable shape, the Tully-Fisher 
relation of field spirals is well fit, the cluster mass-to- light ra- 
tio has the right size, a morphology density relation results, 
the B — V color distribution appears to be well consistent 
with observations, we find luminosity segregation for early 
type galaxies, and we can reproduce the Faber- Jackson rela- 
tion for cluster and field ellipticals. Given the approximate 
treatment of key physical processes, we think that these are 
remarkable successes. However, it is important to note that 
changes in some of our model assumptions can have a strong 
effect on the results. This makes it possible to isolate the con- 
sequences of specific physical assumptions, thereby guiding 
further physical modeling. At this point it is worth noting 
that the number of free parameters in our models is actu- 
ally quite small, in practice, about the same as in direct 
hydrodynamical simulations of galaxy formation. 

One powerful strength of semi-analytic models is that 
they provide the full history of galaxy formation. Using the 
present time to normalize the models, they can make predic- 
tions at high redshift, where observational data can be used 
to put additional strong constraints on physical processes. 
We will discuss the evolution of the galaxy population of our 
cluster in a companion paper (in preparation). 

Combining high-resolution N-body simulations with 
semi-analytic galaxy formation methods is currently the 
only way to simulate directly the formation of a rich clus- 
ter and its constitutent galaxies in their proper cosmological 
context. We think that this approach can be fruitfully ex- 
ploited for study of other aspects of galaxy formation in 
hierarchical cosmologies. 
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